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Smith,  R.E.  1992.  Opus,  An  Integrated  Simulation  Model  for  Transport  of  Nonpoint- 
Source  Pollutants  at  the  Field  Scale:  Volume  I,  Documentation.  U.S.  Department  of 
Agriculture,  Agricultural  Research  Service,  ARS-98,  120  pp. 

Opus  (not  an  acronym)  is  a computer  model  for  the  transport  of  material  in  soil  and 
surface  water.  The  model  is  a simulation  tool  for  studying  the  potential  pollution  from 
various  agricultural  management  practices.  It  simulates  water  movement  that  results 
from  rainfall  and  other  weather  inputs  and  that  is  affected  by  soil,  crop,  topography, 
and  many  types  of  management  actions  and  water  use  influencing  the  surface  condi- 
tion. Opus  includes  models  for  the  growth  of  plants;  development  of  cover,  water  use; 
uptake  of  nutrients;  cycling  of  soil  nitrogen,  phosphorus,  and  carbon;  transport  of 
adsorbed  pesticides  and  nutrients;  interaction  of  surface  water  and  soil  water;  runoff; 
and  erosion.  This  tool  allows  the  user  to  choose  between  (a)  detailed  simulation 
involving  data  on  the  time-intensity  pattern  of  rainfall  and  (b)  a more  lumped  ap- 
proach using  either  recorded  daily  rainfall  or  stochastically  generated  rainfall.  To 
illustrate  its  use.  Opus  is  applied  to  a set  of  rather  comprehensive  data  from  a water- 
shed study  in  Watkinsville,  GA. 

KEYWORDS:  runoff,  erosion,  infiltration,  simulation,  pollution, 
hydrology,  soil  water,  transport 

No  warranties,  expressed  or  implied,  are  made  that  the  computer  programs  described 
in  this  publication  are  free  from  errors  or  are  consistent  with  any  standard  of  program- 
ming language,  or  that  the  programs  will  meet  a user’s  requirement  for  any  particular 
application.  The  U.S.  Department  of  Agriculture  disclaims  all  liability  for  direct  or 
consequential  damages  resulting  from  the  use  of  the  techniques  or  programs  docu- 
mented herein. 

Trade  names  are  used  in  this  publication  solely  to  provide  specific  information. 
Mention  of  a trade  name  does  not  constitute  a guarantee  or  warranty  of  the  product  by 
the  U.S.  Department  of  Agriculture  or  an  endorsement  by  the  Department  over  other 
products  not  mentioned. 

While  supplies  last,  single  copies  of  this  publication  may  be  obtained,  at  no  cost,  on 
request  from  USDA-ARS-NPA,  Water  Management  Research  Unit,  AERC,  CSU, 

Fort  Collins,  CO  80523. 

Copies  of  this  publication  may  also  be  purchased  from  the  National  Technical 
Information  Service,  5285  Port  Royal  Road,  Springfield,  VA  22161. 
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LIST  OF  SYMBOLS 


This  list  contains  the  most  significant  or  the  most  used 
symbols  but  is  not  exhaustive.  Locally  defined  symbols  may  not 
be  included  but  do  not  conflict  with  this  list.  A few  of  the 
symbols,  to  retain  accepted  usage,  have  different  meanings  (as 
indicated  here)  in  different  contexts.  These  symbols  are  used 
only  when  the  meanings  differ  sufficiently  to  avoid  ambiguity 
and  when  usages  do  not  overlap.  In  addition,  common 
abbreviations,  such  as  N for  nitrogen  and  C for  carbon,  are 
used  in  the  text.  Units  throughout  the  text  are  abbreviated  as 
follows : 


m = 

meter 

mm  = 

millimeter 

1 

liter 

gm  = 

gram 

kg  = 

kilogram 

ha  = 

hectare 

min  = 

minute 

Svmbol 

Definition 

Unit 

a 

cross-sectional  area 

m2 

A 

surface  area 

m2 

A 

P 

pond  area  at  h = 0 

m2 

A 

u 

surface  area  contributing  at  the 

upstream  point  of  a channel 

m2 

b 

parameter  used  with  local 

B 

definitions 

parameter  used  with  local 

c 

definitions 

General,  locally  defined 

Ch 

parameter  used  with  subscripts 

conversion  coefficient 

ha-mm/1 

C 

concentration  by  weight 

C 

a 

concentration  of  solute  adsorbed 

to  soil  particles 

kg/kg 

C 

P 

concentration  of  phosphorus 

kg/1 

C 

concentration  of  solute  in 

w 

liquid  soil  water 

kg/1 

VI 
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Symbol  Definition 


C 

s 

C 
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C 

z 

d 


E 

P 


E 

ps 

E 

s 


f 


f 

pw 


f 

rs 

f 

rm 

f 

rs 


concentration  of  sediment  in 
surface  water 

equilibrium  concentration  of 
sediment  in  surface  water  at 
current  conditions 
concentration  of  solute  at  time  0 
rate  of  transfer  of  suspended 
particles  between  surface  water 
and  soil  surface 
potential  rill  detachment  rate 
actual  rill  detachment  rate 
net  rate  of  loss  or  entrainment 
diffusivity  of  soil  water 
soil  thermal  diffusivity 
efficiency  of  cultivation  mixing 
parameter  for  efficiency  of  plant 
respiration 

daily  potential  plant- transpirable 
water 

daily  potential  soil  evaporation 
daily  soil  evaporation 
daily  potential  evapotranspiration 
daily  potential  evapotranspiration 
less  intercepted  water 
rate  of  infiltration 
content  of  soil  carbon  by  weight 
fraction  of  clay  in  surface  soil 
NH^  factor  in  decay  of  residue 
fraction  of  pesticide  on  plant 

surface  that  is  subject  to  washoff 
content  of  soil  residue  by  weight 
fraction  of  plant  residue  that  is 
metabolic 

fraction  of  plant  residue  that  is 


kg/1 


kg/m2/min 

kg/m2/min 

kg/m2/min 

kg/m2/min 

mm2/min 

mm2/min 


mm 

mm 

mm 

mm 

mm 

mm/min 


structural 


Symbol  Definition 


Unit 


f fraction  of  sand  in  surface  soil 

s 


W 

F(  ) 


F 

P 


F 

m 


F 

s 


g 

g 

g 

G 


c 

s 


G 

e 

h 

h 

o 


temperature  factor  in  decay  of 
residue 

water-content  factor  in  decay  of 
residue 

objective  function  for  numerical 
minimization 

maximum  fraction  of  surface  area 
shaded  by  plants 

leaf  area  index:  leaf  area  divided 
by  shaded  area 

index  of  potential  maximum  leaf 
area 

total  leaf  area,  considering  all 
plants 

fraction  of  soil-surface  area 
shaded  by  mulch 
fraction  of  total  shaded  soil 
surface  area 

gravitational  acceleration 
capacity  for  transport  of  sediment 
actual  sediment  transport  rate 
parameter  of  capillary  and  water 
deficit  for  infiltration  model 
value  of  G for  surface,  uncrusted 
soil 

value  of  G for  wetting  front, 
including  crust  effect 
depth  or  height  of  water 
Chapter  5:  normal  flow  depth  for 
surface  waterflow 
Chapter  4:  initial  value  of  h 

m 


= t 

o 


m/sec/sec 

kg/m2/min 

kg/m2/min 

mm 

mm 


mm 

m 


m 


viii 


at  t 


Symbol  Definition 


Unit 


h average  depth  of  superelevation 

of  water  table  above  draintiles 


h^  maximum  height  of  superelevation 

of  water  table  above  draintiles 


H 

H 

c 


H 

o 


i 


I 

I 

a 


I 

P 

j 


J 

k 

x 


k 

o 


K 


K 

oc 


K 

g 


K 

s 


K 

u 


total  soil  water  potential,  y>  - z 

integral  capillary  drive 

infiltration  parameter 

daily  net  effective  radiant 

energy  in  terms  of  equivalent 

evaporated  water 

index  subscript:  represents 

spatial  step  in  Chapters  4 and  5 

depth  of  infiltrated  water 

initial  abstraction  in  SCS  Curve 

Number  formula 

I at  ponding  during  rainfall 

index  of  time  increment  in 

numerical  expressions 

jacobian  tri-diagonal  matrix 

decay  coefficient  of  carbon  pool 

x in  equations  of  residue  decay 

dimensionless  kinematic  flow  number: 

= S L/h 
o o 

local  hydraulic  conductivity  of 
soil 

adsorption  ratio 

base  adsorption  ratio  (multiplied 
by  carbon  content  of  soil  to  get 
K for  gravitational  flow  of  soil 
water 

effective  saturated  hydraulic 
conductivity,  K at  ^ = 0 
K for  diffusive  flow  in  soil  water 


mm 

mm 

mm 

mm 

mm 

mm 

mm 

mm 


day  1 


mm/min 

1/kg 


mm/min 

mm/min 


Symbol  Definition 


M 


M 

M 


lv 


m 


M 


pm 


M 


n 


Pi 


u 


length  of  flow  path  of  surface 
water 

length  of  flow  path  through 
channel 


Unit 


m 


m 


L 

u 

length 

factor  in  USLE  equation 

m 

p 

content 

of  pesticide  in  soil 

kg/ha 

m T 
rL 

content 

of  lignin  in  soil  residue 

g/m2 

mrN 

content 

of  nitrogen  in  soil  residue 

g/m2 

mrC 

content 

of  carbon  in  soil  residue 

g/m2 

mrP 

content 

of  phosphorus  in  soil 

residue  g/m2 

number  of  identical  surfaces 

feeding  into  sides  of  a channel 
(1  or  2) 

plant  dry  matter  in  leaves  and  stems  kg/ha 
density  of  mulch  cover  kg/ha 

total  plant  dry  matter  kg/ha 

total  potential  maximum  plant  dry 

matter  kg/ha 

mass  of  soil  layer  kg/ha 

Chapter  5:  Manning  roughness 
coefficient 

Chapter  7 : content  of  plant 

nitrogen  kg/kg 

number  of  elementary  units  of 

runoff  making  up  the  catchment  area 
probability  function  of  a continuous 
random  variable 

probability  of  a discrete  random 
variable 

total  depth  of  daily  or  storm 

rainfall  mm 

factor  for  cropping  practice  in  USLE 
equation 


X 


Symbol  Definition 


Unit 


s 

c 

S 

S 

o 

S 

u 

t 

T 


flux  of  local  surface  water,  per 

unit  width  m2/min 

flux  of  sink  or  source  in  soil 

column  mm/min 

flux  of  soil  water  due  to  gravity  mm/min 

inflow  to  computational  layer  mm/min 

flux  of  soil  water  into  saturated 

zone  that  feeds  draintiles  mm/min 

outflow  from  computational  layer  mm/min 


peak  rate  of  runoff 
flux  of  suspended  sediment,  by 
volume 

flux  of  soil  water  due  to 
capillary  gradient 
surface  discharge 
inflow  to  pond 
outflow  from  pond 

outflow  of  sediment  from  watershed 
rate  of  rainfall 
hydraulic  radius,  a /p 
vector  of  R values  at  nodes  along 
flow  path 

daily  incoming  solar  radiation 
energy  factor  in  USLE  equation 
parameter  of  soil-water  storage,  used 
in  SCS  Curve  Number  method  for 
runoff 

slope  of  soil-water  characteristic, 

dd /dip 

slope  of  flowing  water  surface 
slope  of  soil  surface  along  flow  path 
slope  factor  in  USLE  equation 
time  (units  vary  with  process) 
temperature 


mm/min 

m3/min 

mm/min 

m3/min 

m3/min 

m3/min 

kg/m2 

mm/min 

m 

m 

langley 

ha-mm/N/hr 


mm 


mm  1 


°C 


Symbol 

Definition 

Unit 

T 

K 

Kelvin  temperature 

°K 

T . 
mm 

minimum  daily  temperature 

°C 

T 

max 

maximum  daily  temperature 

°C 

u 

velocity  of  water 

m/min 

u* 

shear  velocity 

m/min 

U. 

1 

size-3  vector  representing  cross 

correlations  in  weather-generation 

model 

V 

local  variable  of  integration  or 

local  coefficient 

V 

s 

velocity  of  settling  of  particles 

m/min 

V 

volume  of  water 

m3 

w 

width  of  section  of  surface  flow 

(Ch.  5) 

m 

w 

s 

water  equivalence  of  a snowpack  unit 

mm 

W 

width  of  distributed- flow  unit 

(Ch.  2) 

m 

X 

distance  along  path  of  surface  flow 

m 

size  of  particle 

mm 

ys 

spacing  of  subsurface  parallel  drains 

m 

Y 

random  variate , value  0 to  1 

z 

depth  from  surface 

mm 

z 

c 

side  slope  (ratio  of  h to  v)  of 

furrow  or  channel 

zt 

total  depth  of  soil  profile  modeled 

mm 

Z 

r 

depth  of  rooting 

mm 

a 

Chapter  2:  curvature  parameter  in 

soil-water  characteristic  function 

Chapter  3:  shape  parameter  for  gamma 

probability  distribution  function 

P 

Chapter  3:  scale  parameter  for  gamma 

probability  distribution  function; 
elsewhere:  relaxation  fraction  in 


numerical  computations 


Symbol  Definition 


Unit 


7 

7s 

7w 

S(  ) 


A 


e 


6 

e 


e 


K 

S 


\ 


V 


K 

p 

pb 


o 


T 

c 


specific  weight  of  water 
heat  capacity  of  a soil  particle 
volumetric  heat  capacity  of  water 
iterative  correction  term  for 
equation  of  surface  runoff 
slope  of  curve  for  saturation  vapor 
pressure  at  mean  air  temperature 
exponent  in  relation  of  K to  0: 

= (2  + 3A)/A 

content  of  volumetric  soil  water 
content  of  residual  soil  water,  a 

parameter  in  the  soil  characteristic 
function 

content  of  natural  saturated  water, 

6 at  rp  = 0 

normalized  water  content:  = 

(d-6  )/(9  -6  ) 
r s r 

thermal  conductivity  of  surface  soil 
pore-size  distribution  parameter, 
log  slope  of  soil  hydraulic 
characteristic 

rate  coefficient  for  kinetic  model 
of  adsorption  partitioning  exchange 
albedo 

wetted  perimeter 

wetted  perimeter  along  the  bottom 
of  a furrow  or  channel 
specific  gravity  of  a particle 
wetted  perimeter  of  the  sides  of  a 
furrow  or  channel 
standard  deviation  of  a random 
variable 

critical  shear  at  bed  on  which 
erosion  is  assumed  to  occur 


cal/cm3/deg 

cal/cm3/deg 


mcal/cm 


min  1 

m 

m 


m 


N/m2 


Symbol  Definition 


Unit 


rN  stress  factor  for  crop  growth  due  to 

N deficiency  in  soil 

rg  shear  at  bed  of  surface  flow  N/m2 

t stress  factor  for  crop  growth  due 

to  suboptimum  temperature 
stress  factor  for  crop  growth  due 
to  inadequate  soil  water 
<f>  modifying  coefficient  for  crop  effects, 

used  in  various  formulas  to  estimate 
erosion  detachment 

X ^ residual  vector  in  generation  model 

for  daily  weather 

V>  capillary  potential  of  soil  water  mm 

air-entry  parameter  for  soil 

capillary  characteristic  function  mm 
<j>  relative  time  weighting  in  numerical 

method  (0.5  to  1.0)  (0.7  used 
in  Opus) 
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OPUS,  AN  INTEGRATED  SIMULATION  MODEL  FOR  TRANSPORT  OF 
NONPOINT- SOURCE  POLLUTANTS  AT  THE  FIELD  SCALE: 

VOLUME  I,  DOCUMENTATION 

Roger  E.  Smith 
1.  INTRODUCTION 


Models  allow  scientists  to  mathematically  state  what  they  have 
learned  about  nature  and  natural  processes.  The  familiar  equa- 
tions F = ma  and  e *=  me2  are  models.  The  advent  of  computers 
■has  given  scientists  the  ability  to  represent  knowledge  of  the 
interactions  of  several  related  processes.  To  the  extent  that 
the  resulting  model  duplicates  nature,  the  model  allows  the 
study  of  "what  if"  questions  and  the  completion  of  certain 
limited  experiments  that  may  require  years  or  centuries  by 
physical  experiment. 

Several  models  have  recently  been  produced  that  deal  with  the 
many  interactive  processes  occurring  in  the  water-based  ecosys- 
tem of  an  agricultural  field  (SPAW,  NTRM,  CREAMS,  and  EPIC,  to 
name  a few) . All  such  models  are  imperfect  in  some  sense 
because  of  the  limited  knowledge  of  the  various  processes  and 
their  interactions.  Although  good  theories  may  exist  for  parts 
of  this  ecosystem,  at  least  under  laboratory  conditions,  natu- 
ral heterogeneities  always  confound  the  application  of  these 
theories  to  field  conditions.  Other  processes  may  be  described 
only  imperfectly,  and  empiricisms  or  conceptual  approximations 
may  be  required  to  represent  some  parts  of  the  system. 

Nevertheless,  the  need  to  study  the  movement  of  potential 
pollutants  from  today's  chemically  intensive  agriculture  has 
stimulated  the  building  of  and  the  use  of  models  for  agricul- 
tural hydrology  and  transport.  In  all  cases,  published  models 
are  a compromise  between  scientific  rigor  and  the  practical 
limits  as  to  the  kinds  of  data  that  can  reasonably  be  expected 
to  be  obtained  for  any  model  applications. 

Opus  (not  an  acronym)  is  offered  herein  as  a potential  improve- 
ment over  currently  available  models  for  the  processes  of 
agricultural  hydrology.  It  too  is  imperfect  and  has  the  same 
limitations  as  described  above.  However,  Opus  does  overcome 
many  of  the  weaknesses  identified  in  several  existing  models 
(Smith  and  Ferreira  1989),  and  it  has  relatively  advanced 
approaches  to  the  simulation  of  water  movement  and  chemical 
transport . 


R.E.  Smith,  U.S.  Department  of  Agriculture,  Agricultural 
Research  Service,  Water  Management  Research  Unit,  Fort 
Collins,  CO. 
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The  purpose  of  Opus  is  to  allow  the  study  of  the  effects  of 
management  and  weather  inputs  on  the  movement  of  water  and 
potential  pollutants  within  and  from  a small  catchment.  Using 
the  available  physical  and  chemical  information  and  the  guide- 
lines in  the  User  Manual  (vol . II),  one  can  describe  the  plant 
profile,  soil  profile,  and  topography  of  a small  catchment; 
subject  that  information  to  either  historical  or  statistically 
generated  inputs  related  to  weather  and  to  management;  and 
select  from  a variety  of  output  types  and  details  to  summarize 
the  results  of  interest. 

This  volume  I describes  the  theory  and  methodology  used  in  the 
various  processes  simulated  by  the  model,  to  provide  the  user  a 
background  for  understanding  the  capabilities  and  limits  of  the 
model.  A second  volume,  the  User  Manual,  describes  the 
required  input  data  and  the  output  options,  and  gives  (where 
possible)  a guideline  for  the  selection  of  appropriate  parame- 
ters for  a given  condition. 

In  addition  to  and  in  conjunction  with  water  movement,  Opus 
simulates  the  processes  of  sediment  transport  and  chemical 
transport,  and  also  the  cycles  of  carbon,  nitrogen  and  phospho- 
rus involved  in  microbiological  decay  in  soil,  flow  of  heat  in 
soil,  and  growth  of  crops.  The  model  shows  the  changes  in 
these  responses  due  to  a wide  variety  of  simulated  management 
practices.  Simulations  of  runoff,  plant  growth,  and  movement 
of  soil  water  are  always  included;  at  the  user's  option,  any  or 
all  of  the  other  processes  may  be  simulated  concurrently. 

The  components  of  Opus  were  designed  or  selected  to  be  of 
generally  comparable  precision.  In  other  words,  it  was 
intended  that  one  component  not  be  treated  with  a level  of 
detail  incongruent  with  another  component  of  equal  importance 
to  the  simulation  outcome.  Nevertheless,  the  model  components 
involve  various  degrees  of  development  across  the  scientific 
disciplines  represented.  For  example,  the  modeling  of  plant 
growth  and  response  to  stress  has  received  less  overall  atten- 
tion than  has  the  modeling  of  surface  runoff.  There  are  neces- 
sarily some  empiricisms  in  parts  of  the  model;  no  model  is 
entirely  physically  based,  in  the  pure  sense  of  the  term. 

Some  parts  or  processes  in  Opus  have  been  intentionally 
limited,  as  a choice  between  the  complexity  and  the  signifi- 
cance to  most  envisioned  applications.  For  example,  much  of 
the  thermodynamics  of  the  soil - freezing  and  snowmelt  processes 
involves  latent  heat  of  the  ice/water  system,  which  is  not 
considered  in  Opus  except  in  the  case  of  snowmelt.  The  use  of 
degree-day  methods  for  estimating  snowmelt  is  understood  to  be 
very  approximate  for  a small  area  such  as  that  treated  by  Opus. 
Also,  diurnal  variations  in  the  surface -soil  temperature  are 
not  treated:  the  average  daily  surface  soil  temperature  is 


considered  acceptable  for  study  of  long-term  variations  in 
profiles  of  soil  temperature.  This  eliminates  detail  that  may 
make  a limited  improvement  in  long-term  simulation  and  also 
eliminates  the  need  for  great  detail  in  input  data. 

The  scope  of  Opus  is  limited  to  runoff  areas  that  are  charac- 
terized by  a single  soil  profile,  have  hydrologic  complexity 
limited  to  first-order  or  second-order  channel  networks,  and 
contain  a single  cropping  system  or  plant-cover  regime.  This 
does  not  mean  that  only  one  crop  per  year  is  considered  but 
rather  that  any  given  crop  combination  must  represent  the 
entire  field  area  of  interest.  The  model  is  not  limited  to 
annual  crops , and  it  can  treat  perennials  such  as  grassed  and 
grazed  pasture.  Opus  can  simulate  periods  of  a few  days  (with 
an  intensive  look  at  one  or  two  hydrographs)  to  many  decades. 

The  model  allows  several  choices  for  simulation  of  the  hydro- 
logic  process.  The  choice  will  generally  reflect  the  amount  of 
detail  of  data  that  a user  is  expected  to  have.  The  hydrologic 
options  are  diagrammed  in  figure  1.  Implied  in  the  variety  of 
choices  here  is  that  the  model  is  intended  to  have  potential 
uses  ranging  from  management  to  research.  The  available  data 
should  be  appropriate  for  the  intended  use  of  the  model.  Sto- 
chastically generated  weather  may  be  adequate  for  long-term 
simulations  to  observe  alternate  management  schemes,  but  to 
look  at  details  of  hydrologic  processes  such  as  peak  runoff  and 
sensitivity  of  erosion  to  microtopographic  changes,  one  should 
have  data  on  the  pattern  of  rainfall  intensity.  The  physically 
based  simulation  of  infiltration  from  rainfall  is  impossible 
without  such  rainfall  data.  Because  these  data  are  often 
recorded  by  noting  the  times  and  depths  when  the  rate  of  accu- 
mulation changes,  or  the  graph  slope  "breaks,"  these  data  are 
commonly  and  herein  referred  to  as  breakpoint  rainfall  data. 

Much  of  the  model  is  original,  but  parts  have  taken  advantage 
of  process  models  developed  by  others.  Included  is  a daily 
weather  generation  model,  which  was  taken  with  minor  adaptation 
from  the  WGEN  model  of  Richardson  (1981).  The  evaporation  and 
evapotranspiration  potential  is  based  on  the  model  of  Ritchie 
(1972),  which  is  also  used  in  the  CREAMS  and  EPIC  models.  The 
daily  runoff  model  is  basically  a modified  SCS  Curve  Number 
approach  that  was  developed  by  Williams  (Smith  and  Williams 
1980)  and  is  similar  to  that  in  the  EPIC  and  SWWRB  models 
(Williams  et  al . 1984).  The  soil  erosion  model  accompanying 
the  daily  model  is  from  MUSLE  (Williams  and  Berndt  1977). 

The  soil  microbial  system,  including  cycles  of  carbon,  nitro- 
gen, and  phosphorus,  is  a daily  time  step  version  of  the  Centu- 
ry model  (Parton  et  al . 1988a, b).  That  model  was  modified  by 
Parton  for  Opus  to  operate  on  daily  time  scales  and  to  include 
the  addition  of  fertilizers  and  plant  residue  from  the  plowing 
under  of  surface  material. 
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Opus  operates  on  a basic  daily  cycle  of  simulations  (illustrat- 
ed in  fig.  2),  with  an  accounting  cycle  of  1 year  for  purposes 
of  summaries.  In  the  more  intensive  hydrology  option  based  on 
rainfall  intensities,  a smaller  interstorm  cycle  is  used  for 
the  redistribution  of  soil  water  between  storms  when  more  than 
one  storm  occurs  in  a day.  The  processes  of  management 
changes,  plant  growth,  snowmelt,  and  nutrient  cycling  are 
implemented  on  a daily  time  step. 

The  remaining  chapters  in  this  volume  (2  through  9)  detail  the 
methodology  used  in  the  simulated  processes  and  indicate  how 
the  processes  interact.  Volume  II,  the  User's  Manual,  further 
describes  the  structure  and  various  options  of  the  program  and 
gives  information  on  specific  data  requirements.  References 
cited  in  volume  I are  listed  after  chapter  9. 


Figure  1. 

Options  for  method  of  hydrologic  simulation  available  within 
the  Opus  program. 


Read  beginning  and 
ending  dates. 


Stop 


Figure  2. 

Basic  daily  and  annual  cycles  of  continuous  simulation  for 
Opus. 


2.  CHARACTERIZING  THE  PHYSICAL  SYSTEM 


Specifying 

Catchment 

Topography 


Opus  deals  with  the  space  over  and  within  small  catchments, 
from  above  the  plant  to  below  the  root  zone  in  the  soil.  This 
is  almost  never  a simple  system,  and  any  computer  simulation 
begins  with  abstractions  to  describe  the  system  in  terms  that 
allow  satisfactory  simulation  of  the  processes.  Two  of  the 
most  important  geophysical  features  of  a catchment  are  the 
surface  topology  (slopes  and  flow  paths  of  surface  water)  and 
the  soil  profile  (soil  hydraulic  properties  and  how  they  may 
vary  with  depth  from  the  surface) . This  chapter  indicates  how 
these  two  key  hydrologic  components  are  described  for  Opus. 


the  Topography  is  a key  factor  in  the  transport  of  potential  pol- 

lutants within  and  from  a catchment.  Topography  controls  the 
paths  of  surface  waterflows  and  influences  the  velocity  of 
those  that  occur  from  rainfall  runoff  or  from  irrigation. 
Finally,  topography  can  enhance  or  inhibit  the  erosion  of 
surface  material. 

Actual  catchments  are  almost  always  complex  in  shape  at  any 
level  of  detail.  To  simulate  the  surface  hydrology  without 
requiring  excessive  amounts  of  input  and/or  computational 
detail,  it  is  necessary  to  describe  the  catchment  in  geometri- 
cally simple  terms  while  retaining  the  most  important  features 
in  the  hydraulics  of  the  surface  flow. 

The  daily  hydrology  option  is  a lumped- time  approximation  and 
is  necessarily  also  a spatially  lumped  method.  Topographic 
information  is  useful  in  that  method  only  in  general  measures, 
such  as  mean  flow  lengths  and  slopes.  Some  of  the  following 
discussion  will  therefore  not  apply  because  of  the  lumping 
involved.  On  the  other  hand,  the  infiltration-based  hydrology 
option  is  relatively  sensitive  to  topography,  and  the  runoff  is 
routed  over  the  space  of  the  simulated  areas  according  to  the 
specified  geometry  (described  below) . 

The  method  used  in  Opus  to  interpret  the  actual  topography  of  a 
catchment  is  illustrated  in  figure  3.  From  our  present  under- 
standing of  surface-water  hydrology,  we  recognize  the  following 
features  that  should  be  preserved  in  the  geometrically  simpli- 
fied description:  Catchment  area;  mean  path  length  for  shallow 
distributed  flows;  mean  path  length  for  channelized  flows, 
where  applicable;  and  net  slope  of  the  mean  flow  path. 

As  shown  in  figure  3,  the  properties  of  each  part  of  the  flow 
path  are  conserved  separately.  In  addition,  the  mean-slope 
profiles  may  be  preserved  and  treated  in  the  distributed 
hydrology  option.  As  explained  in  the  User  Manual,  the  profile 
is  described  by  specifying  the  slope  at  each  point  along  the 
profile  where  there  is  a change  in  the  rate  of  change  of  slope 
(illustrated  in  fig.  4).  Convergence  of  the  surface  (distrib- 
uted) flow,  as  when  water  flows  toward  one  point  at  the  catch- 
ment boundary,  may  also  be  simulated.  As  illustrated  in  figure 
3,  this  convergence  is  simulated  by  preserving  the  width  at  the 
catchment  divide  in  the  transformed  geometry. 
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Abstracting  the  Real  Catchment 

The  basic  topology  of  the  catchment  in  the  model  is  limited  to 
(a)  a single  distributed  flow  surface  contributing  uniformly 
along  one  side  of  a first-order  channel,  (b)  a matched  pair  of 
such  surfaces  contributing  to  both  sides  of  a first-order  chan- 
nel, or  (c)  item  b plus  a small  area  contributing  to  the  upper- 
most point  of  the  channel.  The  total  simulated  area  may  be 
made  up  of  N units  of  such  topology,  integrated  by  a second 
channel.  The  second  channel  does  not  directly  receive  input  of 
surface  water  from  a runoff  surface,  but  takes  the  outflow  from 
the  N catchments  uniformly  along  its  length.  Such  a channel 
receiving  multiple  inputs  would  simulate,  for  example,  the 
function  of  a terrace  outlet  channel. 

The  structure  described  above  can  be  used  in  varying  forms  and 
degrees  of  approximation  to  represent  most  topologies  found  in 
small  catchments.  Asymmetrical  cases  (where  a first-order 
channel  receives  contributions  from  different- sized  areas  on 
its  two  sides,  or  where  the  slope  profiles  on  each  side  may  be 
significantly  different)  can  be  treated  in  Opus  only  by  taking 
mean  properties  (forcing  symmetry)  in  each  case.  Inclusion  of 
detail  of  simulation  and  necessary  data  to  treat  asymmetrical 
cases  explicitly  is  believed  to  be  uneconomic  for  the  purpose 
of  the  model,  because  of  the  additional  computation  time 
required.  When  the  area  on  one  side  of  a channel  is  quite 
small  compared  with  the  other,  the  small  side  can  (for  most 
purposes)  be  added  to  the  large  side  without  significant  error. 

An  important  feature  to  be  considered  in  describing  the  flow- 
path  topology  for  hydrologic  simulation  is  the  role  of  agricul- 
tural row  furrows  in  controlling  the  path  of  water  across  a 
given  field.  This  path  seldom  corresponds  to  the  path  that 
would  be  taken  by  surface  water  in  the  absence  of  furrows.  In 
the  case  of  contour  or  near-contour  plowing  with  or  without 
terraces,  the  flow  path  may  be  entirely  different.  In  some 
cases  where  management  changes  sufficiently  during  a rotation 
period,  the  flow-path  topology  may  change  from  year  to  year,  or 
may  change  at  the  end  of  a season  if  the  rows  are  obliterated. 
For  this  reason,  the  Opus  model  provides  for  alternate  descrip- 
tions (sets  of  descriptive  variables)  of  the  flow  topology 
within  the  same  field  area.  Case  1 describes  the  flow  topology 
for  flow  unaffected  by  furrows,  as  on  a meadow.  Case  2 
describes  flow-path  geometry  when  the  flow  path  is  affected  by 
a furrow  direction  superimposed  on  the  natural  contours.  If 
the  model  is  used  for  natural  catchments  (for  example,  untilled 
pasture),  only  the  case-1  topology  needs  to  be  described. 

Figures  5 and  6 show  topology  of  two  possible  cases  of  flow 
path  for  one  hypothetical  field.  Subscript  1 is  for  natural 
flow  (case  1),  and  subscript  2 is  for  furrow- directed  flow 
(case  2).  In  figure  5,  flow  follows  the  natural  shape,  flowing 
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across  the  contours  and  thus  converging  toward  the  outlet  (here 
pictured  as  a small  pond).  In  the  methodology  of  Opus,  the 
width  of  the  flow  surface  at  the  upper  end  of  the  flow  path  is 
calculated  from  the  user's  estimate  of  the  mean  length  of  the 
flow  path,  the  field  area,  and  the  width  of  the  flow  path  at 
the  bottom.  If  the  same  field  is  tilled  as  in  figure  6,  case-2 
hydrologic  topography  is  produced,  and  the  furrows  cause  runoff 
to  flow  from  left  to  right  along  the  furrow  "roicrochannels . " 
Convergence  of  flow  cannot  occur  as  in  case  1 (unless  the  fur- 
rows are  overtopped  and/or  eroded  away).  In  the  furrowed  case, 
the  total  flow  path  is  lengthened  and  the  path- slope  profiles 
are  correspondingly  changed.  The  total  change  in  elevation 
from  top  to  bottom  is  the  same,  and  the  topological  parameters 
should  preserve  it.  This  may  be  stated  mathematically,  using 
subscripts  for  the  respective  cases,  as  follows: 

f 1S1(x)dx  -=  f 2S»(x)dx  [1] 

J0  J0 


where  Lx  and  L2  are  total  path  length  for  cases  1 and  2, 
respectively,  Sj  and  S2  are  case-1  and  -2  slopes,  and 
x is  distance  along  flow  path. 

Figure  7 gives  two  other  examples  in  which  simplified  hydrolog- 
ic geometry  is  obtained  from  natural  topography.  In  the  upper 
example,  because  natural  flow  occurs  normally  with  the  con- 
tours, the  field's  actual  shape  should  be  distorted  somewhat  to 
preserve  the  actual  flow-path  length  in  the  simplification. 

This  is  because  the  flow,  as  shown,  is  at  some  angle  to  the 
(constructed)  field  border  and  is  thus  somewhat  longer.  Thus 
the  simplified  width  is  made  slightly  less  than  the  topographic 
width  to  preserve  the  total  area  and  mean  path  length.  Another 
way  to  describe  this  is  to  say  that  the  area  is  represented  by 
a parallelogram  of  the  given  area. 


In  the  lower  case  in  figure  7,  flow  diverges  as  it  moves  from 
the  hilltop  toward  the  channels.  This  case  could  be  simplified 
into  two  planes,  each  contributing  to  the  side  of  a different 
channel,  with  the  two  channels  then  joining.  To  preserve  the 
area  ABC  and  a selected  mean  flow  length  (L  ..  ) , the  width  of 
each  abstracted  plane  is  calculated  as  Wx  =^area/L  The  mean 
width  in  this  case  would  be  less  than  L Opus  accounts  for 

this  situation  by  apportionment  of  surface  flow  along  the 
receiving  channel  length.  On  the  other  hand,  if  the  width  so 
calculated  is  greater  than  L then  converging  flow  is  indi- 
cated, which  is  explicitly  modeled  in  the  equations  for  surface 
routing  (see  ch.  5). 


Terrace  systems  are  common;  figure  8 is  a representation  of  a 

terrace  system.  The  Opus  topological  approximation  requires 

that  the  individual  areas  of  terrace  drainage  be  represented  as 

equal.  Usually  this  is  more  or  less  true.  Further,  it  is 

assumed  that  terrace  spacings  are  identical.  This  is  even  more 

commonly  the  case.  Note  the  two  sequential  channels  labeled 

L . and  L 0 . 
cl  c2 

This  case  illustrates  another  important  aspect  of  hydrologic 
topology.  If  the  rows  are  plowed  closely  parallel  to  the  con- 
tours and  to  the  terraces,  the  major  direction  of  runoff  will 
be  across  the  furrows  and  will  occur  only  after  a certain 
amount  of  surface  water  has  accumulated.  Opus  recognizes  this 
case  by  investigating  the  furrow  slopes.  When  slopes  are  neg- 
ligible and  runoff  exceeds  an  estimated  furrow  storage  based  on 
furrow  depth  and  width,  Opus  simulates  furrow  overtopping. 


Implementing  Topography  in  the  Model 

Once  the  subject  topology  has  been  conceptually  simplified,  it 
must  be  described  to  the  computer  program.  The  minimum 
requirement  for  input  to  Opus  is  description  of  case-1  (unfur- 
rowed) topology.  This  is  sufficient  for  natural  catchments  or 
unplowed  pastures  and  meadows.  Case-2  topology  is  also 
described  whenever  (a)  the  initial  furrow  depth  is  significant 
or  (b)  management  operations  creating  significant  furrow  depths 
are  included  in  the  operations  menu  (described  below).  Undis- 
turbed fields  where  no  management  operations  are  used  need  to 
have  only  the  case-1  topology  described. 

As  noted  above,  from  one  to  three  elements  make  up  the  topolog- 
ical network,  through  which  water  is  routed  in  order.  The 
first  element  is  a flow  surface  across  which  flow  Is  assumed  to 
be  uniformly  distributed,  but  the  surface  is  not  necessarily  a 
plane  surface.  This  element  may  represent  one  of  two  relative 
scales.  At  the  smallest  scale,  where  furrows  or  ridges  are  a 
significant  part  of  the  surface  area  (such  as  large  furrow 
terraces  used  on  severely  sloping  land) , the  first  element  may 
represent  the  sides  of  the  ridges  and  may  be  only  a few  meters 
long.  In  this  case,  the  furrow  or  swale  between  ridges  becomes 
the  channel  element  2.  Where  ridges  are  close  as  in  ordinary 
furrows,  it  is  not  recommended  to  represent  fields  in  this 
detail,  because  it  is  unnecessary  to  route  water  over  a flow 
path  1 meter  or  less  long. 

For  ordinary  use,  it  is  preferable  to  use  a larger  scale 
approach  and  consider  all  the  furrows  on  a slope  together  as 
the  first  element.  The  input  data  for  Opus  allow  the  user  to 
describe  the  flow  geometry  of  the  furrows  without  needing  to 
treat  each  one  as  the  second  (channel)  element  in  the  routing 
sequence . 
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Characterizing 
Soil  Profile 


The  elementary  catchment  unit  of  an  Opus  field  is  made  up  of 
either  one  or  two  (M  - 2)  distributed  flow  elements  contribut- 
ing to  a channel  (element  2).  The  topology  of  Opus,  for  sim- 
plicity alone,  is  characterized  by  symmetry.  Thus  there  can  be 
a number  (N)  of  elementary  catchments,  but  they  must  be  equal. 
An  additional  upstream  area  (A  ) that  contributes  to  the  head 
of  the  channel  (element  2)  may  also  be  designated.  Figure  9 is 
a scheme  of  an  example  system  for  a case  where  M = 2,  there  are 
three  elements  in  the  field  (one  surface  and  two  channels,  not 
counting  the  A areas) , and  four  units  make  up  the  field  (N  = 
4).  U 

If  specified,  the  third  element  (the  second  channel)  receives 
uniformly  distributed  input  from  elements  2 along  its  length. 

It  can  also  have  an  upstream  area  (A  ) contributing  separately 
from  the  N basic  catchment  units.  At  the  third  element  level, 
however,  no  multiples  (N)  make  up  the  total  field;  there  can  be 
only  one  channel  element  3.  Figure  10  shows  the  third  element 
acting  as  a terrace  channel  outlet,  such  as  in  figure  8.  In 
figure  10,  all  A^  values  are  0 and  M = 1. 

Within  the  program,  the  flow  length  of  each  element  is  divided 
into  computational  increments.  There  can  be  from  3 to  20  such 
increments  for  any  element  flow  path,  depending  on  length  of 
the  flow  path.  The  number  of  segments  is  generally  taken  as 
72l;  with  L in  meters. 

Convergence  of  flow  along  element  1 is  found  from  the  geometry 
parameters.  As  discussed  above,  the  mean  flow  length  (L  ) and 
the  area  of  distributed  flow  surface  are  preferentially  pre- 
served, and  the  mean  width  is  found  by  division:  W = A /L  . 

If  the  channel  (element  2)  length  is  greater  than  W,  thEn  ^flow 
is  assumed  parallel  in  element  1.  But  if  the  channel  length  is 
shorter,  a convergence  rate  is  calculated  and  used  in  the  sur- 
face-water routing.  Convergence  is  not  allowed  when  parallel 
flow  along  furrows  of  case  2 has  been  otherwise  indicated. 

More  details  and  examples  of  how  to  specify  various  conditions 
within  this  system  are  given  in  the  User  Manual. 


the  The  soil  profile  is  assumed  to  consist  of  one  or  more  soil 

horizons,  with  each  horizon  differing  from  its  neighbor  in  some 
significant  manner.  Such  differentiation  may  consist  of  (a) 
differences  in  texture  and  thus  hydraulic  characteristics  or 
(b)  significant  differences  in  nutrient  or  other  organic  char- 
acteristics . 

An  Opus  user  may  describe  a soil  profile  made  up  of  as  many  as 
10  horizons.  For  computational  purposes,  the  specified  profile 
is  subdivided  systematically  (by  the  program)  into  smaller  lay- 
ers, as  illustrated  in  figure  11.  Smaller  divisions  are  main- 
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tained  at  shallower  depths,  where  the  most  rapid  soil -water 
changes  occur.  Opus  also  internally  subdivides  the  soil  pro- 
file into  three  sectors  for  the  dynamics  of  residue  and  nutri- 
ents: a litter  layer  at  the  surface,  a microbiologically 

active  soil  layer,  and  the  deepest  layer.  Presumably,  these 
divisions  separate  microbiologically  dissimilar  environments. 
The  active  layer  is  assumed  to  extend  from  the  soil  surface  to 
the  horizon  interface  nearest  a depth  of  200  mm  (Parton  et  al . 
1988b),  or  to  that  depth  if  only  one  horizon  exists. 

The  soil  profile  is  an  important  hydrologic  element  because  it 
is  a porous  medium  through  which  water  flows  and  is  also  a 
medium  in  which  plants  grow  and  extract  water.  Thus  the  soil 
horizons  must  be  described  in  quantitative  terms  that  relate  to 
the  soil's  hydraulic  characteristics.  The  statement  that  a 
certain  horizon  is  a "Bl"  horizon,  for  example,  does  not  give 
information  that  can  be  used  in  simulating  its  water-holding 
and  water-movement  characteristics.  The  same  is  true  to  a cer- 
tain extent  for  a soil's  taxonomic  name  or  color;  that  is, 
those  labels  do  not  indicate  a soil's  hydraulic 
characteristics . 

Whenever  the  soil  is  partially  saturated,  liquid  soil  water 
occurs  in  intergranular  interfaces.  Because  of  the  strong 
surface  tension  of  water,  a significant  negative  pressure, 
which  is  locally  continuous  through  the  soil,  is  characteristic 
of  the  soil  water  at  a particular  water  content.  The  lower  the 
water  content,  the  smaller  the  effective  radii  of  the  water 
surface  at  the  granular  interfaces,  and  the  higher  the  negative 
pressure.  This  pressure  is  usually  measured  in  terms  of  effec- 
tive potential  or  head  of  water.  It  is  the  gradient  of  this 
potential  that  causes  water  to  move  in  the  soil.  Most  manage- 
ment models  today,  by  contrast,  treat  soil -water  movement  using 
concepts  of  storage  filling  and  draining. 

The  function  used  in  Opus  to  describe  the  relation  of  water 
content  to  matric  potential  is  illustrated  in  figure  12.  The 
general  expression,  which  can  describe  the  characteristic 
relation,  may  be  stated  as 

e - [i  + W\)a)'x/a  [2] 

in  which  9 is  normalized  volumetric  water  content,  defined  as 

e = (0  - 0r)/(6s  - 6r)  [3] 

and  tj)  is  matric  soil-water  potential  (mm)  , 

is  a parameter  representing  air-entry  potential  (mm)  , 
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a is  a curvature  coefficient  affecting  the  shape 
of  the  curve  near  A, 

A is  a pore-size  distribution  parameter, 

9 is  volumetric  water  content, 

6 is  residual  water,  theoretically  the 
r water  unextractable  by  capillary  forces,  and 

8 is  water  content  at  0 matric  tension, 
s 

This  expression  is  similar  to  that  of  van  Genuchten  (1980)  but 
relates  directly  to  the  Brooks  and  Corey  relation  (1964),  which 
contains  parameters  that  represent  physical  concepts. 

The  companion  relationship  describing  the  unsaturated  hydraulic 
properties  of  soils  is  the  relation  of  hydraulic  conductivity 
(K)  to  water  content.  One  robust  expression  for  this  relation 
is 


K K 0£  [4] 

s 1 J 

in  which  K is  effective  saturated  (if>  - 0)  hydraulic 
conductivity  (L/T) , 

and  according  to  Brooks  and  Corey  (1964) , e may  be  approximated 
as 

e - (2  + 3A)/A  [5] 

In  general,  is  smaller  for  soils  with  large  particle  sizes 
and  is  larger  for  soils  with  fine  particle  sizes.  Similarly,  A 
is  small  for  soils  with  a well-distributed  range  of  particle 
sizes  and  is  large  for  soils  with  a very  uniform  particle  size, 
such  as  sand.  The  value  of  K varies  widely  but  follows  a 
trend  similar  to  that  of  A. 

To  describe  a soil  horizon,  therefore,  Opus  requires  five  major 
parameters:  A,  0g,  9 , anc*  K ■ In  Opus,  a is  fixed  for 

simplicity  at  5.0?  T^e  selectionSof  reasonable  values  for  the 
hydraulic  characteristics  of  a particular  soil  horizon  is  an 
important  task  for  the  user  of  a physically  based  model.  Opus 
takes  advantage  of  the  extensive  data  assembly  and  analysis  by 
Rawls  et  al . (1982)  to  give  the  user  a survey  of  the  expected 
values  of  these  critical  parameters.  This  use  of  regression 
provides  a way  to  estimate  reasonable  parameters  in  most  cases 
from  the  user's  knowledge  of  the  soil  texture  (Brakensiek  et 
al . 1984).  The  accuracy  of  Opus  in  simulation  depends  in 

large  measure  on  the  accuracy  of  the  parameters  chosen.  Exper- 
imentally determined  soil  data  are,  of  course,  most  preferable. 
Other  sources  of  such  soil  hydraulic  data  exist  (e.g.,  Holtan 
et  al . 1968).  The  importance  of  characteristic  relations  in 
the  simulation  of  soil -water  movement  is  discussed  in  chapter 
4. 
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Figure  3. 

The  geometric  variables  used  to  describe  topography  of  a 
catchment.  The  descriptive  variables  represent  a geometric 
simplification  of  the  real  catchment,  which  retains  the  major 
hydrologically  sensitive  measures  of  the  catchment.  These 
include  area,  mean  flow  path  of  surface  water  (L  ),  channel 
length  (L  ),  and  convergence.  p 


Figure  4. 

Varying  slope  profile  of  a catchment  described  by  specifying 
each  point  at  which  slope  trend  changes,  assuming  each 
slope  segment  to  be  an  arc.  For  example,  the  slope  between 
x2  and  x„  changes  smoothly  and  needs  no  intermediate 
specification  points. 
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Figure  5. 

Example  geometry  of  a convergent  field  with  flow  not 
controlled  by  furrows.  A mean  length  is  estimated  as  shown. 
The  upper  catchment  width  would  be  very  large  compared  to 
the  width  at  bottom,  to  represent  convergence  of  flow  lines 
toward  the  pond. 


Figure  6. 

The  same  field  as  in  figure  5,  plowed  to  produce  very 
different  flow  path.  The  width  and  length  are  both  changed, 
but  the  area  remains  the  same.  Now  a channel  of  length  Lc| 
receives  water  from  the  furrows. 


Figure  7. 

Two  more  examples  of  how  surface  flow  is  specified  for 
fields  of  different  topological  characteristics.  Divergent  flow 
from  the  lower  field  is  simulated  by  two  different  approximate 
planes.  L is  average  flow  length. 


Figure  8. 

Representation  of  a typical  terrace  system.  Terrace 
systems  are  common.  Opus  provides  for  a second-order 
channel  (length  L-)  to  collect  flow  from  the  terrace 
intercept  channels  (length  Lc|). 
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Figure  9. 

Example  scheme  of  field  topography.  Opus  limits  complex 
geometries  to  symmetrical  patterns,  as  shown  here.  Opus 
allows  upslope  areas  that  feed  into  the  upper  end  of  chan- 
nels, here  labeled  A „ for  first-order  channels  and  A . for 
second-order  channels.  The  overall  area  is  divided  info  four 
units,  and  each  unit  feeds  from  both  sides  of  the  channel 
(M  = 2). 


Ap 
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Figure  10. 

The  terrace  channel  system  as  described  for  Opus,  in  simple 
geometry.  Each  of  the  eight  units  feeds  to  one  side  of  its 
receiving  channel  (M  = 1). 


Figure  1 1. 

Soil  horizons  (e.g.,A,  B,  C)  are  specified  by  the  user  and 
further  divided  into  computational  increments  within  the 
program. 
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CAPILLARY  POTENTIAL,  4',  mm 

Figure  12. 

The  unsaturated  soil  hydraulic  characteristics  are  repre- 
sented by  an  analytic  expression  that  allows  fitting  to  most 
soil  types.  Water  content  is  normalized  to  relative  saturation, 
O,  as  shown  in  equation  [3], 
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3.  METEOROLOGICAL  PROCESSES 


Alternatives  for 
Weather  Simulation 


The  field  unit  simulated  by  Opus  may  be  thought  of  as  an  eco- 
system, including  the  plants  and  soil  down  through  the  root 
zone.  In  this  conceptual  framework,  the  meteorological  pro- 
cesses are  the  driving  forces  for  this  system,  and  the  descrip- 
tion of  those  processes  is  a crucial  part  of  the  simulation 
process.  The  meteorological  variables  that  "drive"  or  excite 
the  system  here  include  daily  solar  radiation  (measured  in 
langleys) , daily  maximum  temperature,  and  precipitation.  These 
"drivers"  may  be  input  or  in  some  cases  may  be  generated  from 
input  statistical  parameters  representing  an  analysis  of  local 
historical  records. 


Whenever  available,  records  of  daily  maximum  and  minimum  tem- 
peratures, daily  radiation,  and  rainfall  for  a location  are 
input  to  drive  the  simulation.  The  input  format  is  described 
in  the  User  Manual  (vol . II).  In  practice,  all  required  data 
(on  meteorological  variables)  may  not  be  available  for  many 
locations  for  which  simulation  studies  are  desired,  or  may  not 
have  been  collected  frequently  enough  to  allow  physically  sound 
simulations.  Most  likely  the  user  will  have  only  Weather 
Bureau  summary  data  perhaps  including  monthly  means  of  tempera- 
tures, or  will  have  only  daily  precipitation  records,  or  will 
have  some  combination  of  these  data.  So  Opus  has  been  designed 
to  accept  a variety  of  combinations  of  such  actual  records, 
including  daily  values  or  historical  monthly  means  of  any  of 
the  maximum  temperatures,  minimum  temperatures,  and  net 
radiation . 

The  weather  model  in  Opus  is  built  around  a variation  of  the 
WGEN  model  (Richardson  and  Wright  1984) , which  requires  as  a 
minimum  no  additional  data,  but  falls  back  on  a set  of  fur- 
nished parameter  maps  for  the  continental  United  States.  In 
addition,  Opus  accepts  data  on  monthly  average  temperatures  and 
radiation,  and  reconstructs  either  a smoothly  changing  record 
or  a serially  correlated  random  record  with  the  given  charac- 
teristics. Thus  a variety  of  available  data  can  be  used  in 
conjunction  with  the  WGEN  model  to  produce  the  necessary 
weather  record  for  simulation. 

Fundamental  to  this  model  is  the  idea  that  these  weather  param- 
eters are  intrinsically  linked- -that  temperature  variations  and 
radiation  variations  from  day  to  day  should  be  related  to  each 
other  and  to  the  occurrence  pattern  of  rainfall.  The  model  can 
maintain  both  the  statistical  features  of  observed  time  series 
and  the  generally  observed  interdependencies  between  physically 
related  variables. 

The  data- generating  process  and  its  assumptions  are  summarized 
below.  In  many  cases,  although  daily  temperatures  and  radia- 
tion values  are  available  for  a certain  period  of  record,  the 
user  may  choose  to  generate  data  when  longer  simulations  are 
desired,  instead  of  creating  and  handling  massive  files  of 
input  data. 
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Simulation  of  a 

Precipitation 

Record 


When  either  daily  or  pattern  ("breakpoint")  data  for  rainfall 
are  available  for  the  period  for  which  the  simulation  is  to  be 
performed,  the  data  are  given  to  the  model  in  the  form  of  a 
sequential  data  file.  If  such  data  are  not  available  or  if  it 
is  desired  to  simulate  for  periods  longer  than  the  available 
record,  a statistical  model  of  the  occurrence  and  depth  of 
daily  rainfall  may  be  used  to  generate  rainfall  sequences 
appropriate  to  the  local  climate.  Note  that  only  daily  rain- 
fall can  presently  be  stochastically  generated;  breakpoint  data 
must  be  provided. 


Simulating  Occurrence  of  Daily  Rain 

The  occurrence  of  daily  rainfall  is  modeled  by  WGEN  as  a Markov 
chain  process,  simulating  the  binary  sequence  of  wet  days  and 
dry  days.  On  rainfall  days,  or  wet  days,  the  amount  is 
randomly  distributed  according  to  a two-parameter  gamma 
distribution.  A wet  day  is  defined  as  one  on  which  more  than 
0.25  mm  (0.1  inch)  of  rain  falls.  The  following  description  is 
a summary  of  the  WGEN  methodology,  and  the  reader  is  referred 
to  Richardson  (1981)  or  Richardson  and  Wright  (1984)  for  a more 
complete  discussion. 

The  Markov  chain  model  of  a binary  process  deals  with  transi- 
tion probabilities,  that  is,  the  probability  of  the  transition 
from  a wet  day  yesterday,  for  example,  to  a dry  day  today.  Let 
p,(w|d)  represent  the  probability  of  a dry  day  following  a wet 
day  on  day  i.  Then  clearly, 

Pi(w|w)  = 1 - pi(d|w)  [6] 

and 

pi(d|d)  = 1 - pi(w|d)  [7] 

Thus,  specifying  two  of  these  four  transition  probabilities, 
p(d|w)  and  p(d|d),  fully  defines  the  model.  These  probabili- 
ties change  with  seasonal  climatic  changes.  With  enough 
historical  data,  one  can  get  a good  picture  of  the  change  in 
transition  probabilities  in  some  climates  over  periods  as  short 
as  5 days  to  1 week.  For  purposes  of  Opus  (and  in  WGEN)  it  is 
considered  adequate  to  define  weather  trends  by  month.  Thus  i 
is  used  in  equations  [6]  and  [7]  as  a monthly  index,  and  the 
rainfall-occurrence  model  requires  a total  of  24  parameters. 
Richardson  and  Wright  (1984)  have  tabulated  these  parameters 
for  many  locations  in  the  United  States.  These  data  help  in 
modeling  for  ungaged  basins  and  are  reproduced  in  the  User 
Manual . 
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Simulating  Depth  of  Daily  Rainfall 

The  probability  of  the  amount  of  rain,  v,  on  wet  days,  p(v) , is 
estimated  by  a two-parameter  function  for  gamma  probability 
density: 


a-1  -v//3 

p(v)  = — — , v,  a,  (3  > 0 [8] 

r (a) 

The  symbols  a and  fi  are  the  parameters  for  shape  and  scale, 
respectively.  Like  the  transition  probabilities,  they  are 
assumed  to  vary  by  month  at  each  location.  Stochastic  simula- 
tion of  a daily  amount  involves  generating,  on  days  for  which 
runoff  is  predicted  to  occur,  a random  variate  Y:  0<Y<1.  Then 
Y is  made  to  have  the  distribution 


Y(v)  = V p(v)  dv  [9] 

J0 

which  is  inverted  to  solve  for  the  precipitation  value  v satis- 
fying this  relation.  When  equations  [8]  and  [9]  estimate  val- 
ues of  v less  than  0.25  mm,  they  are  ignored  and  another  value 
is  produced.  This  procedure  is  consistent  with  the  definition 
of  a wet  day  given  above. 


Simulation  of  Daily  The  simulation  of  daily  radiation  (in  langleys ) , daily  maximum 
Radiation  and  temperature,  and  daily  minimum  temperature  depends  on  options 

Temperature  chosen  by  the  user  and  on  the  existence  of  any  input  records  of 

those  variables.  The  basic  model  for  the  variation  of  these 
three  variables  (heretofore  referred  to  as  "weather"  variables) 
assumes  that  in  the  mean  they  follow  a sinusoidal  annual  varia- 
tion, have  a random  variation  about  this  mean  pattern,  are 
cross-correlated,  and  are  correlated  with  the  wet/dry  pattern 
of  rainfall.  This  part  of  the  model  is  also  taken  from  the 
WGEN  model  of  Richardson  (1981)  and  is  briefly  described  below. 
Richardson's  methods  have  been  slightly  modified  to  allow  oper- 
ation in  the  Celsius  scale  of  temperatures.  The  changes  to 
this  basic  weather  model  to  account  for  recorded  input  data  are 
described  in  the  following  section. 


Generating  Mean  Values  for  Daily  Temperature  and  Radiation 

If  monthly  data  on  temperature  and  radiation  are  provided,  the 
simplest  method  for  obtaining  estimates  of  record- averaged 
daily  weather  values  is  to  interpolate  between  the  monthly  mean 
values  obtained  from  the  input  data,  assuming  that  the  monthly 
mean  applies  on  the  day  at  the  middle  of  the  month.  This  pro- 
vides a smoothly  changing  value  for  each  record  of  maximum  and 
minimum  daily  temperatures  and  daily  solar  radiation.  With 
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this  method,  the  value  for  these  parameters  on  a given  date 
will  be  identical  for  each  year.  For  some  simulation  purposes, 
this  estimate  of  weather  data  will  be  adequate.  For  this  and 
other  options,  Opus  applies  a bias  upward  or  downward,  depend- 
ing on  whether  the  day  is  dry  or  wet,  respectively. 


Generating  a Stochastic  Record  of  Daily  Temperature  and 
Radiation 

To  obtain  a more  realistic  sequence  of  daily  weather  variables, 
with  or  without  dependence  on  recorded  monthly  mean  values,  an 
autoregressive  generation  scheme  may  be  used  for  each  of  these 
weather  variables.  Monthly  mean  values  may  be  used,  but  daily 
values  are  generated  based  on  the  model  described  by  equations 
[6]  through  [9]  above.  Parameters  are  obtained  from  Richard- 
son's parameter  maps  of  standard  deviations  and  variable 
cross -correlations  (in  app . D of  the  User  Manual).  Thus  daily 
radiation,  for  example,  can  be  a random  variable  that  has  the 
mean  and  standard  deviation  as  specified  by  the  input  data,  but 
also  is  correlated  with  yesterday's  value.  The  value  for  daily 
radiation  shifts  downward  if  rainfall  occurs. 

The  procedure  for  generating  daily  values  of  T , T . , and  R. 
(maximum  temperature,  minimum  temperature,  and ' solar " radiation , 
respectively)  was  presented  by  Richardson  (1981)  . Each  of 
these  daily  values  is  assumed  to  be  related  to  its  value  yes- 
terday in  an  autoregressive  manner  as 

xt(j)  = Axi_1(j)  + Bx-(j)  [10] 


where 


x. ( i ) is  the  vector  of  residuals  of  T (i  = 1) , T . 
l,.  , . max  mm 

(j  = 2),  and  Ri  (j  = 3)  on  day  i, 

X . is  a size-3  vector  of  independent  random  components, 

1 and 

A and  B are  3X3  matrices  whose  elements  define  the  serial 
and  cross-correlation  properties  of  these  variables. 

Richardson  (1981)  found  that  although  means  and  standard  devia- 
tions vary  in  space  and  time,  it  is  a very  acceptable  approxi- 
mation to  assume  constant  values  of  correlation  relations 
described  by  matrices  A and  B.  Opus  uses  such  fixed  matrices. 
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Since  x.  is  a residual  value,  with  mean  0 and  standard  devia- 
tion 1.6,  the  generated  value  of  any  given  variable,  v^j),  is 
then  found  as 

vi ( j ) = v.(j)  + xi(j)ai(j)  [11] 

where  j indices  are  as  defined  above, 

v.  is  the  mean  value  of  v for  day  i,  and 

l 

o.  is  the  standard  deviation  of  v for  day  i. 

The  generation  model  distinguishes  two  sets  of  parameters  for 
each  location:  one  is  for  wet  days  and  one  is  for  dry  days. 
Each  set  includes  parameters  for  the  mean  and  amplitude  of  the 
Fourier  two- term  series  description  of  the  mean  value  of  each 
variable,  plus  a mean  and  amplitude  for  a Fourier  model  of  the 
coefficient  of  variation  of  each.  The  mean  and  amplitude  of 
the  coefficient  of  variation  for  radiation  are  assumed  the  same 
for  all  locations,  based  on  analysis  of  the  U.S.  data.  In 
operation,  Opus  converts  coefficient  of  variation  into  standard 
deviation,  to  avoid  possible  division  by  zero  when  operating  in 
the  Celsius  temperature  scale. 


Using  Recorded  Monthly  Data 

When  monthly  mean  recorded  data  are  available  for  any  of  the 
three  weather  variables,  Opus  adjusts  generated  values  to  have 
monthly  mean  values  equal  to  the  corresponding  recorded  means. 
The  model  uses  recorded  data  in  comparison  with  the  Fourier 
means  for  each  month,  and  determines  a shift  factor  that  keeps 
the  generated  mean  values  consistent  with  the  input  values  of 
monthly  means.  The  weather-generating  parameters  may  be  taken 
from  parameter  maps  referred  to  above,  and  are  supplied  as 
input  to  Opus  in  any  case. 


Reading  Recorded  Daily  Data 

When  recorded  daily  data  are  available,  an  option  in  Opus 
allows  the  user  to  bypass  the  weather- generating  model  and  to 
read  actual  recorded  weather  data  on  a day-by-day  basis.  A 
separate  input  file  is  used  for  this  option.  Such  recorded 
weather  data  may  include  the  three  weather  variables  (described 
above),  the  runoff,  and  the  daily  sediment  production.  This 
option  is  included  for  users  wishing  to  simulate  only  ancillary 
processes  such  as  plant  growth,  transport  of  nutrients,  or 
transport  and  decay  of  pesticides. 


Potential 

Evapotranspiration 


The  values  of  temperature  and  radiation  on  a given  day  can  be 
used  to  calculate  the  potential  evapotranspiration  (E  ) on  that 
day.  Actual  transpiration  plus  soil  evaporation  is  less  than 
or  equal  to  E , and  depends  on  plant  and  soil  conditions.  The 
method  used  in  Opus  to  calculate  E is  a modified  form  of  the 
equation  used  in  CREAMS.  It  was  developed  by  Ritchie  (1972), 
starting  from  the  relations  of  Penman.  The  equation  and  its 
components  are 

E - (1  + c ) AH  /(A  + 0.68)  [12] 

two 


where 


c is  a coefficient  expressing  effects  of  wind  and 
humidity , 

A is  slope  of  curve  for  saturation  vapor  pressure  at 
mean  air  temperature,  represented  as 

A = 5304.0  exp  (21.255-5304. 0/T,.)/T  2 [13] 

1\  K. 

in  which 

T is  daily  mean  Kelvin  temperature 
K. 

The  value  of  c is  0.28  in  Ritchie's  equation,  but  Opus  allows 
the  user  to  modify  it  to  account  for  local  conditions.  Con- 
ceptually, humid  conditions  reduce  this  factor,  and  windy 
conditions  increase  it  (Ritchie  et  al . 1976). 

The  second  variable  in  the  E equation  is  defined  as 

Hq  = R.(l  - 0/58.3  [14] 


where 


£ is  albedo  of  field  surface  (described  below) , and 

R^  is  incoming  solar  radiation,  in  langleys/day 

In  these  expressions,  R.  is  the  daily  solar  radiation  in 
langleys/day.  Equation  [14]  is  an  estimating  algorithm  for  the 
daily  net  effective  radiant  energy  (H  ) . The  units  of  E are 
millimeters  of  water,  and  E^_  can  be  thought  of  as  a derived 
daily  variable  for  weather.  Note  that  1 langley/day  = 0.041868 
mj/m2/day . 


21 


The  field  albedo  (f  ) is  calculated  as  a combination  of  area- 
weighted  factors  for  plant  and  mulch  cover,  snow  cover,  and 
surface  soil.  The  albedo  of  dry  surface  soil  is  an  input 
parameter.  Snow  albedo  is  presumed  to  vary  with  time  (in  days) 
after  snowfall  (t  ),  as  follows: 


'snow 


.4  + .5 


e<-V12) 


[15 


Soil  albedo  varies  with  normalized  moisture  content  (0)  and 
compares  with  the  dry  soil  albedo  (£^  ) as  follows: 


’ soil 


W1-0 


0.4  02) 


Plant  and  mulch  albedo  is  assumed  constant  at  0.23. 


[16] 


Time  Distribution  of  PE 

Figure  13  schematically  illustrates  how  Opus  assumes  the  major 
weather  variables  to  be  distributed  in  time.  Time  distribution 
is  important  because  computational  time  intervals,  which  vary 
with  the  hydrologic  conditions,  may  consist  of  any  portion  of  a 
day.  On  days  that  include  a rainfall  event,  the  beginning  and 
ending  times  of  that  event  are  found  or  assumed,  and  the  E for 
that  day  is  assumed  to  be  that  occurring  in  the  nonrainfall 
portion  of  the  daylight  hours.  Radiation  is  assumed  to  be  dis- 
tributed sinusoidally  between  dawn  and  dusk  (as  shown  in  fig. 
13)  on  rainless  days.  From  equations  [12]  and  [14],  this 
implies  a sinusoidal  distribution  of  E as  well.  A small  time 
interval  near  dusk  or  dawn  contributes  considerably  less  incre- 
mental than  does  an  equal  time  near  noon. 

Division  of  PE  Between  Plant  and  Soil  Surface 

For  an  estimated  daily  E^,  the  soil  and  plant  conditions  deter- 
mine the  daily  proportion  expected  to  be  evaporated  from  the 
soil  surface  or  transpired  through  the  plant  leaves,  or  in 
winter  the  amount  evaporated  from  snow  surfaces.  In  addition, 
some  of  the  E^  will  be  used  when  water  is  intercepted  from 
precipitation  and  stored  on  plants  or  surface  material.  Such 
intercepted  water  has  first  access  to  E^..  The  remaining  E^, 
which  we  here  refer  to  as  E ',  is  divided  among  soil  and 
plants . 


Transpiration  is  estimated  on  the  basis  of  total  area  of 
leaves.  That  portion  of  E^'  consumed  by  plants  is  estimated  as 


E 

P 


where 


E 'F  c 
t p s 


[17] 


Fp  is  relative  surface  area  shaded  by  plants,  and 

c is  relative  effectiveness  of  plant  cover  as  a function 
of  leaf  area  index,  as  follows: 


cs  - 1/  [1  + (2 . 5/Fl) s ] o • 2 [18] 

where  F is  leaf  area  index. 

L 

Leaf  area  index  is  defined  as  the  ratio  of  the  total  leaf  sur- 
face of  a plant  to  the  ground  area  projection  or  shade  under 
the  plant.  Equation  [17],  through  F , also  accounts  for  sparse 
plant  cover  (when  plants  do  not  cove?  the  entire  ground  sur- 
face). The  value  of  c can  vary  from  0 to  1 . Further  discus- 
sion of  the  use  of  plant  water  is  in  chapter  7. 


Soil  evaporation  is  estimated  from  the  portion  of  soil  surface 
not  shaded  by  plants  or  surface  litter  (mulch) . The  total 
shaded  area  includes  both  live  plant  cover  and  surface  mulch 
cover,  and  is  estimated  as 


F - F + (1  - 0.5F  )F 
s p p m 

where  F is  relative  area  covered  by  mulch, 
m J 


[19] 


The  relative  amount  of  E ' assigned  to  soil  evaporation  is 
controlled  by  the  soil-water  content  in  the  near-surface 
region.  Without  this  latter  control,  the  potential  soil 
evaporation  part  of  E^'  (E  g)  is 


Eps  - VexP(-l-2Fs>  [20] 

Soil  evaporation  proceeds  in  the  two- stage  process  assumed  by 
Ritchie  (1972)  and  used  in  CREAMS  (USDA  1980).  The  first  stage 
is  short  and  proceeds  at  the  potential  rate  (E  ) based  on 
micrometeorologic  conditions.  The  second  stage^represents 
soil-limited  evaporation,  and  the  rate  is  reduced  proportional- 
ly to  the  square  root  of  time,  provided  soil  water  can  move 
toward  the  surface  at  this  rate.  This  process  reduces  the 
actual  soil  evaporation  (E  ) below  the  potential  represented  by 

E . S 

ps 
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From  equations  [18]  through  [20],  E + E do  not  always  neces- 
sarily sum  to  E . If  the  sum  is  grlater^  the  deficit  is  allo- 
cated so  that  tfie  sum  becomes  E . Of  course,  the  plant  may 
also  not  be  able  to  transpire  at  its  potential  due  to  limited 
soil  water  in  its  root  zone,  but  this  is  not  reflected  in  any 

increase  in  the  soil-surface  evaporation.  If  E + E is  less 

* s p 

than  or  equal  to  E , no  adjustment  is  required. 


Figure  13. 

Time  distribution  of  major  weather  variables  in  Opus. 
Radiation  is  assumed  to  vary  with  solar  angle  as  shown,  and 
temperature  is  assumed  to  follow  a daily  sinusoidal  cycle. 
Radiation  is  reduced  during  a rain,  and  rain  also  reduces  the 
daily  maximum  temperatures. 


4.  FLOW  AND  TRANSPORT  BY  SUBSURFACE  WATER 


Simulating  Flow  of 
Unsaturated  Water 


Soil-water  dynamics  involves  most  of  the  precipitation  falling 
on  the  agricultural  catchment  and  is  thus  a key  element  in  the 
water  transport  of  pollutants,  both  directly  and  indirectly. 

The  zone  of  unsaturated  flow  extends  from  the  soil  surface  down 
to  the  water  table.  In  most  cases,  the  solution  region  in  Opus 
will  extend  to  only  a few  centimeters  below  the  maximum  root 
depth.  In  the  case  of  high  water  tables  with  artificial  drain- 
age, the  solution  region  will  extend  to  the  lower  boundary 
control,  which  is  usually  the  drain  tile  depth. 


Opus  uses  Richards'  equation  to  describe  flow  in  the  unsaturat- 
ed soil  profile.  This  mathematical  expression  is  based  on  a 
dynamic,  one -dimensional  form  of  soil-water  conservation  with 
Darcy  flux,  expressed  as 


a± 

at  az 


where 


[21] 


1 9 is  volumetric  soil  water  content, 

t is  time  (min) , 

z is  depth  from  surface  (mm) , 

q is  flux  (mm/min) , and 

q is  local  inflow  (min  1 ) . 

e 

The  value  q is  assumed  to  be  described  by  Darcy's  law,  which 
relates  flux  to  gradient  of  total  potential  (H)  as  follows: 

q - -K(0)“  ’ [22] 


where 


K is  hydraulic  conductivity  (mm/min) , and 
H is  total  soil  (gravitational  and  matric)  potential 
(V>  + z)  so  that 


mi  = d±  . 

3z  3z 


[23] 


Richards'  equation  assumes  that  the  unsaturated,  one -dimension- 
al vertical  flow  of  water  is  unaffected  by  the  counterflow  of 
air.  The  equation  does  not  attempt  to  account  for  other  soil- 
water  -movement  processes  involving  thermodynamics,  such  as 
vapor  movement  or  local  soil  freezing.  These  assumptions  are 
not  considered  restrictive  in  the  context  of  the  principal 
applications  envisioned  for  Opus,  and  in  comparison  with  other 
sources  of  error. 
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Simulating  Redistribution  of  Soil  Water 


All  common  numerical  solutions  to  equation  [21]  divide  the  soil 
depth  (z)  into  finite  increments  within  which  V and  therefore  9 
and  K(0)  are  considered  locally  uniform.  One  of  the  recurring 
problems  in  the  numerical  solution  of  equation  [21]  is  charac- 
terizing an  effective  local  value  of  K(0)  for  flow  between  the 
adjacent  finite  increments  making  up  the  numerical  grid.  K (9) 
is  a highly  nonlinear  function  and  is  often  changing  rapidly  in 
space  and  time,  whether  or  not  the  soil  texture  changes  between 
layers.  Haverkamp  and  Vauclin  (1979)  have  summarized  and  com- 
pared several  published  techniques  for  K(0)  averaging  within 
numerical  solutions. 


For  Opus,  with  relatively  large  z increments,  the  problem  is 
even  more  acute,  and  the  method  used  takes  unique  advantage  of 
the  analytic  knowledge  of  unsaturated  gravitational  soil  flow. 
Taking  equations  [22]  and  [23]  together,  we  can  consider  the 
flux  (q)  as  being  made  up  of  two  flux  components  corresponding 
to  the  two  components  of  the  total  potential  H (diffusive  flux, 
q^,  and  gravitational  flux,  q ) as  follows: 


q 


q + q 

u Hg 


■k  (e)j^ 

u dz 


V° 


[24] 


The  two  different  subscripts  on  K(0)  (u  = diffusive  and  g = 
gravitational)  represent  two  different  methods  of  calculating 
the  effective  K between  finite  difference  nodes.  The  gravita- 
tional q is  always  downward  (positive  z) . From  kinematic 
analysis®(Smith  1983b),  K for  flow  between  nodes  i and  i+1  is 
best  approximated  for  numlrical  solution  by  heavily  weighting 
toward  K(0)^,  the  K at  the  upper  node. 


For  diffusive  flux,  an  efficient  value  for  K (i,i+l)  can  be 
found  by  assuming  that  ip  varies  uniformly  between  z^  and  z^+^. 
Thus,  while  dip/dz  is  uniform  between  nodes  and  is  continuous 
within  a given  soil  horizon,  dd /dz  will  nevertheless  be  discon- 
tinuous when  soils  change  at  the  layer  interface  between  i and 
i+1.  This  is  illustrated  in  figure  14.  Since  K(0)  is  a 
relatively  simple  relation  (equation  [4]),  the  net  K is  calcu- 
lated as  the  first  moment  of  K(0)  over  the  computational  inter- 
val. To  describe  the  case  of  changing  soils  at  the  interface, 
0,.  is  defined  as  water  content  at  the  j side  of  interface  b, 
with  j = i or  i+1.  Then  the  first  moment  of  K(0)  between  i and 
i+1  is 


K ( i , i+1)  = 
u 


^7  r h 


(v)dv  + 


i b , i J 9 


b , i 


9 

v-f 


b , i+1 


b , i+1  i+1  J 9 . 1 

l+l 


^i+i (v) dv 


[25] 
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With  the  soil's  relation  for  K( 0 ) characterized  simply  as 
described  in  chapter  2,  the  solution  of  equation  [25]  is  rela- 
tively straightforward. 


Using  the  conceptual  division  of  q from  equation  [24],  equation 
[21]  becomes 


s (0) . f^i  - t^[K  ^ + K ] . + 
c l 3t  8 z 1 u 3z  gi 


‘ei 


[26] 


where  s (0)  is  the  local  slope  of  the  hydraulic  characteris- 
tic, i.e.,  d0/d0.  This  method  is  termed  the  Separable  Flux 
(SF)  method  of  numerical  treatment. 

The  first  step  of  the  solution  procedure  for  Opus  is  a time 
step  (At)  selection.  This  is  accomplished  by  calculating  q in 
equation  [24]  for  the  flow  between  all  nodes  at  existing  condi- 
tions, and  choosing  a At  by  examining  the  resulting  net  rate  of 
change  of  0 in  each  computational  layer.  The  critical  condi- 
tions that  are  examined  include  (a)  drying  to  near  0 at  the 
surface  and  (b)  potential  (artificial)  gradient  reversals  that 
would  arise  for  large  time  steps  with  locally  large  differences 
in  q.  An  iterative,  partially  linearized,  implicit  method  is 
used,  which  estimates  values  of  K that  are  fixed  for  the  time 
step,  but  provides  mass -conserving  correspondence  between  the 
new  array  of  values  and  the  new  values  of  0 . This  means  that 
values  of  s (0)  are  chosen  as  chord  slopes  between  old  and  new 
points  on  tfte  0(V>)  relation,  and  are  kept  correct  across  the 
time  step  by  iteration.  The  nonlinear  nature  of  the  0(V>) 
relation  and  its  reversal  of  curvature  require  the  use  of  a 
relaxation  technique,  and  the  relaxation  factor  becomes  larger 
for  cases  in  which  several  iterations  are  required.  Solution 
proceeds  with  the  automatic  selection  of  time  steps  until  the 
time  interval  of  soil-water  redistribution  (not  longer  than  1 
day)  is  exhausted.  Because  of  the  vast  array  of  soil  condi- 
tions met  and  also  the  severe  nonlinearity  of  the  characteris- 
tic relations,  time  steps  can  and  do  vary  as  needed  between 
0.01  and  1440  minutes. 

Local  inflow  (q  ) is  negative  for  both  evaporation  in  the  sur- 
face increment  and  plant  root  extraction  in  all  lower  incre- 
ments. The  total  root-water  use  is  distributed  according  to 
the  water  potential  seen  by  the  plant  at  each  level,  which  is 
proportional  to  z . - ip.  . Continuous  accounting  is  made  of  the 
available  water  in  the  root  zone.  On  a daily  basis,  this  is 
compared  with  the  water  demand  of  the  plant  and  with  the  esti- 
mated "wilting  point"  water  content,  to  limit  plant  uptake  in 
response  to  stress.  (This  is  discussed  in  more  detail  in  ch. 
7.)  Surface-soil  evaporation  is  estimated  by  the  modified 
method  of  Ritchie  (1972)  (ch.  3).  Surface-soil  evaporation  is 
computationally  limited  when  the  surface  nears  15-bar  water 
content,  by  the  soil's  ability  to  furnish  water  to  the  surface 
through  upward  flux. 
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The  lower  boundary  condition  for  equation  [26],  as  indicated 
above,  is  located  either  at  two  computation  nodes  below  the 
root  zone  or  at  the  drain  depth  when  field  drains  are  simula- 
ted. In  the  former  case,  0 is  allowed  to  rise  at  the  boundary 
in  a kinematic  manner  when  a 9 "wave"  passes  below  the  roots 
(positive  d^>,/dz)  . Conversely,  xl>  may  decrease  as  this  "wave" 
passes  out  of  the  root  zone  (negative  d\/>/dz)  . At  the  boundary, 
\J)  is  not  allowed  to  decrease  below  the  value  that  represents  a 
hydrostatic  profile  condition  relative  to  the  water  table. 
Periods  of  relative  drought,  or  large  root- zone  demands,  can 
therefore  cause  upward  movement  from  water  in  the  soil  below 
the  root  zone.  This  potential  flux  is  reduced  with  increasing 
depth  to  the  water  table.  If  draintiles  are  included,  the 
unsaturated- zone  solution  includes  the  saturated  region  above 
the  drain  tiles. 

In  model  simulation,  soil  water  is  redistributed  over  periods 
of  about  1 day.  When  a day  includes  a storm,  the  soil -water 
simulation  period  will  go  to  only  the  beginning  of  the  storm. 
The  next  period  will  include  the  remainder  of  the  day  plus  the 
period  to  6:00  p.m.  the  next  day,  or  to  the  start  of  the  next 
storm,  whichever  comes  first.  The  time  of  6:00  p.m.  is  an 
arbitrary  division  point  for  nonstorm  days.  Equation  [26]  is 
not  solved  during  storm  input;  it  is  not  felt  necessary  or 
appropriate  to  use  this  method  in  lieu  of  the  physically  based 
infiltration  option.  In  any  case,  flux-rate  information  is  not 
available  when  the  daily  hydrology  option  is  used.  Both 
hydrology  options  estimate  the  total  runoff,  total  infiltrated 
depth,  and  time  base  of  the  event. 


Simulating  Rain  Wetting  Profiles 

When  runoff  occurs,  saturation  of  surface  soil  must  have 
occurred,  and  a first-order  estimate  of  the  resulting  wetted- 
soil  profile  is  step  input  of  depth  equal  to  the  infiltrated 
amount  divided  by  (9  -9.).  The  prestorm  value  of  9 is  9.  at 

each  depth.  When  runoff  cloes  not  occur  for  the  daily  option, 
the  effective  average  rate  of  rainfall  is  estimated  by  the 
amount  of  rain  divided  by  the  amount  needed  for  runoff  to 
occur.  This  ratio,  less  than  1,  can  be  used  as  a relative  rate 
(r  ) , and  a profile  water  content  can  be  estimated  from  using 
this  value  of  r^  in  an  inversion  of  the  K(0)  relation  (equation 
[4]).  This  provides  an  estimated  0(r)  for  the  step  profile 
estimate.  It  extends  to  that  depth  necessary  to  provide  for 
the  amount  of  infiltrated  rainwater.  In  either  case  the  esti- 
mated profile  becomes  the  starting  point  for  subsequent  soil- 
water  redistribution.  For  the  breakpoint  hydrology  option, 
when  runoff  does  not  occur,  the  rainfall  data  provide  an  aver- 
age value  for  r,  which  is  used  as  above  to  estimate  a 8 from 
the  inverted  K(0)  relation,  for  the  surface-wetted  profile. 


The  relative  efficiency  of  the  solution  method  in  computation, 
using  the  large  z increments,  is  exemplified  in  figure  15.  The 
initial  condition  profile  is  shown  as  a dashed  line.  The  SF 
method  described  here  is  compared  with  a precise,  fine-mesh 
solution  by  an  implicit  finite-difference  Crank-Nicholson  meth- 
od (Smith  1970).  Given  the  vast  difference  in  finite  increment 
scale,  the  SF  method  is  quite  a faithful  simulation  compared  to 
the  fine-scale  solution,  which  takes  about  2 orders  of  magni- 
tude more  computation  time.  What  is  not  demonstrated  here  is 
the  additional  fact,  with  comparable  large-scale  increments, 
that  the  use  of  common  K-averaging  schemes  and  standard 
numerical  methods  (such  as  the  Crank-Nicholson  scheme)  would 
not  produce  an  acceptable  approximation  to  the  behavior  of 
equation  [ 21 ] . 


Simulating 

Draintile 

Control 


When  draintile  management  methods  are  used  in  areas  of  a high 
water  table,  the  action  and  effect  of  the  draintiles  as  lower 
boundary  controls  on  the  unsaturated  flow  region  may  be  approx- 
imated using  an  analytical  expression.  Figure  16  illustrates 
the  approach  used  to  link  the  saturated  and  unsaturated  flow 
regions.  The  superelevation  of  the  water  table  above  the  par- 
allel drains  is  assumed  to  describe  a half -ellipse  in  shape. 
Thus  the  mean  superelevation  of  the  water  table  (h)  can  be  made 
a simple  function  of  the  maximum  superelevation  (h  ) . From  an 
expression  of  Bouwer  and  van  Schilfgaarde  (1963) , the  flux 
through  the  draintiles  is  a quadratic  function  of  this  maximum 
superelevation  (h  ).  The  unsaturated  flux  from  above,  q(z,t), 
at  level  h , can  be  assumed  constant  over  the  period  At  of 
interest  (usually  equal  to  the  period  of  unsaturated- flow 
redistribution) . It  is  hereafter  termed  q and  can  be  found 
from  the  unsaturated- flow  solution.  The  controlling  water- 
balance  equation  for  the  basic  portion  of  the  system  shown  can 
be  written  as 

H - Mm  • 1271 

in  which 


y is  drain  spacing  (mm) , 

q^  is  drain  discharge  (mm3/mm2/min) , and 

(f>  is  effective  soil  porosity  in  drain  region. 


From  Bouwer  and  van  Schilfegaard  (1963),  we  use  Houghoudt's 
quadratic  equation  for  q^,  which  may  be  written  as 


q 


d 
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[28] 
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Transport 

Heat 


in  which  ca  and  are  geometrically  calculable  parameters  as 
follows : 

c = 4K  /y  2 [29a] 

a s s 

c,  = 8K  z /y  2 [29b] 

b s e ' J s 

The  quantity  z is  a function  of  z^  (the  distance  below  drains 
to  a limiting  interface),  and  is  calculated  internally  by  Opus 
following  graphical  results  from  Bouwer  and  van  Schilfgaarde 
(1963)  . 

Equation  [27]  can  be  formed  into  a lower  boundary  condition  to 
use  in  conjunction  with  the  numerical  solution  of  equation  [24] 
for  the  case  of  draintile  control.  This  boundary  condition  is 
active  whenever  there  is  a positive  value  of  H at  the  drain- 
tiles.  Assuming  the  superelevation  to  have  an  elliptical 
shape,  h is  related  to  H by  a constant  ratio  of  4/tt  . The 
expression  is  partially  linearized  for  solution  in  the  tridiag- 
onal matrix  by  using  h (j-l)h  ( j ) in  place  of  h 2(j).  The 
sumbols  j and  j-1  refer  here  to  time -step  indices  for  current 
and  previous  times,  respectively.  This  allows  the  quadratic 
equation  to  be  put  in  the  same  form  as  the  finite  difference 
equations  for  solution.  Very  little  accuracy  is  lost  in  the 
linearization,  because  in  almost  all  practical  cases,  the 
second-order  term  is  2 orders  of  magnitude  smaller  than  the 
linear  term. 

Flow  toward  active  roots  and  surface  evaporation  can  cause  the 
water  table  for  some  cases  to  drop  below  the  drain  level.  Opus 
treats  this  condition  as  it  does  a profile  without  drains.  The 
relation  between  water  deficit  below  the  drains  and  the  nega- 
tive H at  the  drains  is  assumed  to  be  equal  to  the  relation  of 
a soil  column  in  equilibrium  with  the  potential  H at  the  drain 
level.  This  varies  as  necessary  for  the  solution  of  the  bal- 
ance in  equation  [26].  When  the  water  table  again  rises  to  the 
drain  level  due  to  leaching  from  above,  the  boundary  condition 
reverts  to  the  drain  equation  [27].  Use  of  a fixed  but  very 
small  positive  value  for  s when  ^ > 0 allows  Opus  to  success- 
fully solve  for  positive  as  well  as  negative  matric  heads. 


of  Soil  Simulation  of  the  biological  processes  in  the  soil  of  a field 

area  requires  knowledge  of  the  temperature  at  each  depth  and 
time.  Both  plant  growth  and  microbial  activity  are  temperature 
dependent,  and  these  processes  are  involved  in  several  aspects 
of  the  transport  and  decay  of  pollutants.  Although  large 
changes  occur  through  the  day  in  the  near -surface  region  of  the 
soil  profile  (due  to  diurnal  temperature  changes  in  the  air  and 
diurnal  variation  in  soil  radiant  heating) , Opus  focuses  on  the 
simulation  of  variations  in  mean  daily  temperatures.  Heat 
moves  through  the  soil  by  convective  transport  with  water  and 
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by  diffusion  from  temperature  gradients  induced  by  temperature 
changes  at  the  soil  surface.  In  Opus,  these  changes  are 
treated  in  a manner  similar  to  the  movement  of  water,  with 
thermal  redistribution  between  storms  treated  in  consideration 
of  diffusive- type  flows,  and  rapid  waterborne  heat  convection 
during  a storm  treated  as  if  it  were  a single  pulse  flow. 


Diffusion  of  Soil  Heat 


The  soil  profile  is  a significant  heat-storage  mass,  and 
changes  in  daily  surface  air  temperature  cause  diffusion  of 
heat  into  and  out  of  the  soil.  The  general  one -dimensional 
equation  for  heat  diffusion  is 


3T  _ tw  *^2T 

3t  )st2 


dz- 


[30] 


with  boundary  conditions 

T(0,t)  = T , 0 < t < 1440  min,  and 

T(L, t)  - T°  , 
av 

and  initially 

T(z,0)  = T(z)  arbitrary 


in  which 


is  temperature  (°C) , 

is  maximum  depth  of  simulation, 

is  thermal  diffusivity  (mm2/min) , 

is  time  (min) , 

is  depth  (mm) , 

is  soil  density  (gm/cm3), 

is  volumetric  soil-water  content,  and 

is  average  annual  temperature. 


Soil  thermal  diffusivity  (D^)  is  theoretically  a function  of 
the  mineral  composition,  the  shape  and  orientation  of  soil 
particles,  and  water  content  (de  Vries  1966).  Mineral  thermal 
diffusivity  can  be  computed  approximately  from  knowledge  of 
clay  content,  assuming  randomly  shaped  and  oriented  particles. 
Following  the  procedure  for  this  computation  from  de  Vries 
(1966)  and  assuming  that  the  soil-water  content  is  nonuniform 
but  constant  over  the  time  steps  of  interest,  D^,  can  be  deter- 
mined for  each  soil  computational  layer,  and  the  difference 
form  of  equation  [30]  can  be  quickly  solved  at  each  time  inter- 
val with  a direct  Gauss-Seidel  method.  This  T is  found  at  the 
intersection  of  the  soil  computational  layers,  and  T at  z = L 
is  taken  for  simplicity  as  the  annual  mean  air  temperature.  At 
depths  of  a few  meters,  this  is  a reasonable  assumption  for  the 
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lower  boundary  condition.  Although  equation  [30]  could  be 
solved  through  the  day  to  simulate  the  diurnal  excursions  in 
near-surface  temperatures,  that  detail  is  not  necessary  for  the 
purposes  of  simulating  processes  that  depend  on  soil  tempera- 
ture over  the  time  scale  of  interest  in  Opus. 


Convection  of  Soil  Heat  With  Water 


Heat  may  be  input  to  the  soil  with  the  influx  of  water  at  the 
soil  surface,  and  heat  may  be  transported  within  the  soil  by 
the  movement  of  soil  water.  Complete  thermal  mixing  is  assumed 
within  each  incremental  layer  of  the  flowing  water,  and  initial 
soil-water  temperature  is  assumed  to  be  in  equilibrium  with 
soil  mass.  The  convective  heat  budget  for  a control  volume  of 
soil  may  be  written  as 


7 (T.q. 
'w  ini 


Toqo)  “ dt 


(7  V + 7 V )T  ] 
w w s s oJ 


[31] 


in  which 


T.  - temperature  of  inflow  water  (°C) , 

T = temperature  of  outflow  and  control -element  water 
° (°C) 

7 = heat  capacity  of  water  (cal/cm3/deg) , 

V - volume  of  control-element  water  (cm3), 

7g  = volumetric  heat  capacity  of  soil  (cal/cm3/deg) , 

Vg  - volume  of  soil  in  control  element  (cm3), 
q^  = flux  of  inflow  water  (cm3/min) , and 
qQ  = flux  of  outflow  water  (cm3/min) . 

The  left  side  of  equation  [31]  is  the  net  input  of  heat,  and 
the  right  side  is  the  change  in  heat  in  the  control  element. 


The  general  solution  to  this  linear  differential  equation  may 
be  written  as 


T 
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Ti  + 


(T 


oz 


'V  +bV  +AV\*Vwz/AV 
wz  s | 

V +bV  J 
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[32] 


in  which 


V 

b 

At 

AV 


oz 

wz 


initial  value  of  T , 
initial  water  in  element  at  t 

w 

t - t (min) , and 
(q^  -Zqo)At  = AqAt. 


= t 


For  the  special  case  AV  -*■  0 , one  obtains  the  following  from 
equation  [ 31 ] : 


T = T.  + (T  - T. )exp [ -qAt/(V  + bV  ) ] 
o 1 oz  1 M/  w s J 


[33] 


Equation  [32]  is  solved  for  each  computational  element  in  the 
soil  column  in  which  water  moves  during  a time  step.  Values 
for  the  initial  and  final  water  contents  are  obtained  from  the 
waterflow  calculations.  The  most  dramatic  heat  convection  oc- 
curs in  the  influx  of  generally  cooler  rainwater.  In  the 
absence  of  actual  data,  rainwater  temperature  is  assumed  to  be 
equal  to  the  minimum  daily  temperature  for  the  day  of  occur- 
rence . 


Transport  of 
Solutes 
Through  Soil 


Opus  includes  the  capability  of  calculating  the  transport  of 
partially  soluble  and  partially  adsorbed  chemicals  in  flow  of 
unsaturated  soil  water.  One  of  two  kinds  of  solute  adsorption 
may  be  simulated:  (a)  Adsorption  assumed  to  be  characterized 
by  equilibrium  chemistry.  This  adsorption  is  described  by  a 
linear  adsorption  isotherm.  (b)  Adsorption  assumed  to  be 
kinetic,  with  the  rate  of  transfer  between  adsorbed  and  solu- 
tion phases  a linear  function  of  the  current  phase  ratio  com- 
pared to  the  equilibrium  ratio. 

Equilibrium  adsorption  assumes  an  invariant  equilibrium  ratio 
between  the  solute  adsorbed  on  the  soil  solids  and  the  solute 
in  solution  in  the  soil  water.  Equilibrium  adsorption  can  also 
be  taken  to  mean  that  the  kinetic  transfer  rate  (discussed 
below)  is  very  rapid  compared  with  the  processes  of  soil-water 
movement.  The  equilibrium  ratio  is  termed  K^,  and  is  expressed 
as 

adsorbed  concentration  . 

d G concentration  in  water 
w 


Since  C has  units  of  kg/liter,  and  C has  units  of  kg/kg,  the 
ratio  is  in  liters  per  kilogram.  Further,  the  for  pes- 
ticides is  assumed  to  be  directly  proportional  to  the  local 
content  of  organic  carbon  (f^) , as  follows: 


K 


f,K 
b oc 


[35] 


with  K ^ being  a pesticide-dependent  absorption  constant.  If 
this  formulation  were  universally  accurate,  K would  depend 
only  on  the  type  of  chemical.  Although  some  variation  remains, 
this  expression  does  explain  much  of  the  observed  variation  in 
adsorption  properties  of  many  chemicals.  Values  of  K sug- 
gested by  studies  for  many  pesticides  are  given  in  tfie  User 
Manual  (vol . II). 


The  assumption  of  kinetic  adsorption  accepts  the  concept  of 
but  specifies  a rate  of  transfer  between  the  dissolved  phase 
and  the  adsorbed  phase  as  being  proportional  to  the  extent  to 
which  the  existing  ratio  differs  from  the  equilibrium  ratio. 
Expressing  the  transfer  as  the  change  in  concentration  of 
adsorbed  material  (C  ) per  unit  time,  this  rate  may  be  repre- 
sented in  the  kinetic  assumption  as 
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[36] 


dC 
a 

dt 


i/(C  K, 
w d 


- Ca> 


in  which  u is  rate  coefficient  with  units  of  minutes  1 . 


Transport  With  Equilibrium  Adsorption 

A definition  sketch  for  the  transport  of  solutes  through  an 
arbitrary  numerical  layer  in  Opus  is  given  in  figure  17.  If 
equilibrium  is  assumed  to  always  hold  between  the  adsorbed  and 
dissolved  materials,  then  equation  [36]  reduces  to  C = K^C^, 
and  a differential  mass  - conservation  equation  is  written  for 
the  total  amount  of  material  in  the  layer  as  follows: 


di  Vv  + chKdMs> 


C.Q. 

11 


C Q 
o o 


[37] 


in  which 

C is  concentration  in  solution  (kg/liter) , 
is  adsorption  coefficient  (liter/kg), 

M is  mass  of  soil  layer  (kg/ha) , 

c^  is  a conversion  factor:  ha-mm/liter  (0.0001), 
q is  flow  in  or  out  of  a layer,  and 

C is  concentration  in  or  out  of  a layer,  designated 

by  subscripts  as  follows: 
i is  subscript  representing  inflow,  and 

o is  subscript  representing  outflow. 


With  complete  mixing  in  each  layer  assumed,  and  when  q and  q. 
are  both  positive,  C = C , and  one  term  of  equation  [37] 
disappears.  Expanding  the  differential  on  the  left,  noting 
that  dV/dt  is  AQ,  and  grouping  terms,  equation  [37]  becomes 


dC  q.dt 

w = ni 

C.-C  _ V +Aqt+K,M  c, 
l w wz  n d s h 


[38] 


This  equation  also  holds  for  the  case  of  reversal  of  both  in- 
flow flux  and  outflow  flux,  by  symmetry.  When  both  flows  are 
into  the  layer,  this  equation  will  not  hold,  but  C can  be  found 
by  simple  summation,  because  the  layer  is  a sink.  When  both 
flows  are  out  of  the  layer,  their  concentrations  and  the 
source -layer  concentration  will  be  unchanged  over  the  time 
step . 


34 


Solving  equation  [38]  for  nonzero  Aq , the  equation  for  change 

in  local  solute  concentration  (C)  after  a time  (At) , beginning 

at  time  (t  ) at  which  C = C , is 
z z 


C (At)  - C.  + (C 
1 z 


ci> 


/V  +K  ,M  c,+AV 
I wz  d s h 

l V +K  M c, 

\ wz  d s h 


-V./AV 


[39] 


in  which 

= volume  of  water  entering  layer  (mm) , 

AV  = change  in  water  volume  in  layer  during  At,  and 

C - solute  concentration  at  initial  time  (t  ). 

z z 


Again,  uniform  mixing  of  layer  water  and  inflow  water  is  as- 
sumed within  the  layer.  Values  of  q and  q.  are  obtained  from 
calculation  of  soil-water  movement  for  the  soil  profile  over 
the  period  At.  The  At  for  equation  [39]  is  the  same  time 
interval  as  that  for  soil-water  redistribution  or  storm  inflow, 
and  is  not  greater  than  1 day  (as  described  above) . 


The  outflow  concentration  from  any  layer  may  vary  with  time, 
according  to  equation  [39].  To  simplify  computation  over  long 
periods  of  time,  outflow  concentration  from  a layer  is  averaged 
over  the  computation  interval,  to  become  an  input -flux  concen- 
tration for  the  adjacent  layer.  While  equation  [39]  may  be 
integrated  over  time  to  obtain  the  mean  value,  it  is  more 
straightforward  to  obtain  the  time-average  outflow  concentra- 
tion by  mass  balance  as  follows: 


C 
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C(At)[V(At)+KdMsCh] 


in  which 


[40] 


C — mean  outflow  concentration  over  At, 

V^At)  = liquid  volume  in  layer  at  end  of  At, 

C(At)  — concentration  in  layer  at  end  of  At,  and 
-=  initial  liquid  volume  in  layer 

For  the  special  case  of  through  flow  without  volume  change,  or 
q^  = q , equation  [38]  simplifies  and  the  solution  for  C(At)  is 

C(At)  = Ct  + (Cz  - Ci)exp[-qAt/(V(At)  + K^c^  ] [41] 


The  solution  procedure  begins  at  a layer(s)  that  is  a "source" 
(has  only  outflow)  and  proceeds  from  each  source  in  either 
direction  until  meeting  a "sink"  layer  (which  has  only 
inflows) . 
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Transport  With  Kinetic  Adsorption 


The  option  for  transport  with  kinetic  adsorption  assumes  that 
the  transfer  of  material  between  adsorbed  and  dissolved  phases 
is  not  instantaneous  and  follows  equation  [36].  Because  of 
this  assumption,  a mass-balance  equation  including  both 
adsorbed  and  dissolved  materials  must  be  written.  Figure  17 
illustrates  the  linked  systems.  For  the  dissolved  phase,  one 
can  write 


d_ 

dt 


dC 

(C  V)  = -c,— if  + q.C. 
v w h dt  l 


q C 
no  o 


[42] 


where  V is  assumed  to  be  a steadily  increasing  or  decreasing 
function  of  time  during  the  increment  of  computation 

This  equation  can  be  solved  in  tandem  with  equation  [36]  for 
kinetic  gain  or  loss  of  adsorbed  material.  The  form  of  the 
solution  varies  somewhat  with  the  various  conditions  of  inflow 
and  outflow.  Special  cases  exist  at  the  boundary,  where  evapo- 
rative outflow  will  transport  no  dissolved  material  and  where 
input  may  be  free  of  solutes. 

In  general,  when  the  terms  of  the  equation  are  of  equal  order, 
the  Runge-Kutta  method  can  be  used  in  a piecewise  manner  over 
time  intervals  as  great  as  1 day  or  more  to  simulate  the  kinet- 
ic transfer  of  partially  adsorbed  solutes.  However,  many  pes- 
ticides of  interest  are  highly  adsorbed  to  organic  material  and 
clays  in  the  soil,  and  for  realistic  values  of  u,  the  relative 
change  of  C over  a computational  interval  can  be  considered 
negligible.  In  that  case,  the  equation  reduces  in  general  form 
to 


!hc„v> 


c,  i/(K  .C 
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+ q1c1 
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[43] 


Here  C is  written  for  C , which  is  true  for  complete  layer 
mixing  except  when  q is  evaporative  outflow.  Even  for  slowly 
changing  C , this  equation  can  be  superior  in  stability  to  a 
Runge-Kutta  solution.  The  general  solution  to  this  linear  dif- 
ferential equation,  starting  with  C for  interval  At,  is  writ- 
ten as 


C (At)  = C,  + (C 
w b z 


\b/Aq 


[44] 


where  is  a steady- state  concentration  defined  by 


q.C.  + i/c,  C 
c = i i h a 

b q.  - q + i/c,K 


[45] 


h d 
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and  b in  the  exponent  of  equation  [44]  is  equal  to  the 
denominator  on  the  right  side  of  equation  [45].  For  small 
values  of  Aq  (which  make  the  exponent  in  equation  [44]  very 
large) , the  solution  reduces  to 

Cw(t)  = Cb  + (Cz  - Cb)exp(-bAt/Vz)  [46] 

For  either  equation  [44]  or  [46],  one  can  obtain  the  mean  out- 
flow concentration  for  the  interval  by  integrating  as  was  done 
in  obtaining  equation  [40] . Opus  uses  one  of  these  two  equa- 
tions when  is  larger  than  40,  and  selects  equation  [46]  when 
the  exponent  (b/Aq)  exceeds  10  in  equation  [44] . In  either 
case,  the  time  steps  used  follow  the  solution  appropriately, 
and  the  small  changes  in  C that  occur  with  the  large  values  of 
are  computed  stepwise  by  a mass-balance  equation.  The 
kinetic  option  is  useful  for  only  values  of  v smaller  than 
0.01,  because  larger  values  are,  on  a daily  basis,  practically 
indistinguishable  from  the  equilibrium  case. 


Implementation  by  Computer 

Simulating  the  ensemble  of  many  soil  computational  layers  over 
a time  interval  requires  consideration  of  one  or  more  cases, 
depending  on  direction  of  flow  and  solute  transfer.  A layer 
may  act  as  (a)  a sink,  with  flow  entering  from  both  above  and 
below,  (b)  a source,  with  flow  moving  away  in  both  directions, 
(c)  a flow-through,  with  water  flowing  in  one  side  and  out  the 
other,  (d)  a surface  "source,"  with  solutes  leaving  at  the 
bottom  and  only  water  leaving  at  the  top,  (e)  a surface  flow- 
through, with  water  and  solute  entering  below  but  only  water 
moving  out  at  the  top,  or  (f)  the  reverse  flow  of  case  e.  For 
each  case,  equation  [37]  or  [42]  will  be  modified  by  elimina- 
tion of  terms  or  zero  coefficients.  Solution  of  either  the 
kinetic  equation  or  equilibrium  equation  is  implemented  when 
the  equation  has  been  properly  resolved  for  the  local  condi- 
tions. For  each  case,  the  solution  order  must  be  determined 
properly  in  advance,  proceeding  in  both  directions  from  each 
source  layer.  Information  on  inflow  concentrations  from 
adjacent  layers  is  not  required  for  only  source  layers.  Thus, 
input  values  are  determined  in  advance  for  every  layer.  The 
solution  finds  the  new  local  concentrations  of  both  adsorbed 
chemical  and  dissolved  chemical,  and  solves  for  the  time- 
averaged  value  of  outflow  concentration  for  use  in  the  adjacent 
layer  solution.  The  kinetic  - solution  algorithm  uses  internally 
determined  time  subdivisions  when  dC/dt  is  a rapidly  changing 
function  over  the  interval  of  interest. 
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DRAINTILE  SIMULATION 


Dynamic  balance  of  saturated  zone: 

V -<f>f  hmdy  « ys 

unsaturated  flux  = q(z,t) 
dV  . dh 


Figure  14. 

Diagram  of  changes  in  'I'  and  O across  the  interface  of  two 
different  soils,  and  values  used  at  adjacent  computational 
nodes. 
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Figure  15. 

Comparison  of  the  Separable  Flux  numerical  method  of 
Opus  (with  large  depth  increments)  with  a numerical  solution 
of  the  standard  Richards'  equation  formulation  (with  very 
small  depth  increments).  The  SF  method  does  a good  job  of 
simulating  the  general  distribution  of  soil  water.  Use  of  very 
large  increments  with  standard  numerical  methods  produces 
biased  estimates  of  soil-water  redistribution. 


Figure  16. 

Estimation  of  drainage  for  perched  or  high  water  table  with 
draintiles.  Approximate  equations  taken  from  Bouwer  and 
van  Schilfgaarde  (1963)  are  used  as  a lower  boundary 
condition  for  the  unsaturated  zone  redistribution  calculations. 
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Figure  17. 

Diagram  for  model  of  dissolved/adsorbed  chemical  move- 
ment using  the  soil  increments  and  interlayer  fluxes  from  the 
soil  water  flow  model.  This  figure  shows  the  major  terms 
used  in  the  equations. 
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5.  SURFACE-WATER  FLOW  PROCESSES 


Much  of  the  pollution  from  agricultural  areas  is  carried  in 
runoff  from  precipitation  events.  Surface  runoff  interacts 
with  the  soil  and  transports  dissolved  or  suspended  pollutants 
off  the  field  or  catchment  area.  Although  many  of  the  poten- 
tial pollutants  may  be  within  the  soil  profile,  they  are  often 
concentrated  near  the  soil  surface.  Clearly  a useful  simula- 
tion of  the  pollution  from  a field  requires  an  accurate 
practical  simulation  of  the  movement  of  water  on  the  field 
surface  as  well  as  within  the  soil. 

The  difficulty  in  simulating  surface-water  processes  is  not  so 
much  in  our  knowledge  of  the  water -movement  mechanisms,  but 
rather  in  having  data  appropriate  to  our  ability  to  represent 
the  mathematics  of  these  processes.  Specifically,  even  when 
data  on  topography  and  soil  are  abundant,  knowledge  of  only 
daily  totals  or  even  storm  totals  of  rainfall  is  inadequate  for 
accurate  simulation  of  either  infiltration  rates  or  runoff 
rates.  A storm  of  given  depth  could  produce  large  amounts  of 
runoff  for  high- intensity  rainfall  but  produce  no  runoff  if  the 
same  depth  of  rain  comes  at  low  intensity. 

Nevertheless,  it  is  most  often  the  case  that  situations  in 
which  simulation  is  a useful  tool  do  not  provide  the  modeler 
with  the  rainfall- intensity  data  that  are  necessary  for  using 
infiltration  methods.  For  this  reason,  Opus  (like  CREAMS  in 
this  respect)  allows  the  user  to  choose  a surface -hydrology 
method  that  is  appropriate  to  lumped,  daily  rainfall  data.  In 
addition,  since  the  processes  of  surface-water  transport  are  so 
dependent  on  surface -water- flow  information,  the  daily  hydrolo- 
gy option  also  uses  a "lumped"  method  for  estimating  erosion 
and  sediment  production  from  the  catchment.  The  user's  choice 
of  hydrologic  option  will  depend  on  both  the  quality  of  the 
data  available  and  the  particular  objectives  of  the  model  user. 

Because  of  the  limited  information  contained  in  daily  rain 
summaries  and  also  the  lumped  nature  of  the  resulting  runoff- 
estimation  methodology,  the  information  on  field  topology  (ch. 
2)  finds  only  limited  use  in  the  daily  hydrology  option. 
Nevertheless,  advantage  is  taken  of  all  information  from  the 
field  topology  that  can  be  used  to  improve  the  accuracy  of 
daily  runoff  prediction.  This  improvement  includes  a refined 
estimate  of  the  relative  time  of  concentration,  which  is  used 
to  estimate  the  peak  runoff  rate. 


Simulation  of  Runoff  The  minimum  hydrologic  information  needed  for  Opus  from  a pre- 
From  Daily  Rainfall  cipitation  event  is  the  amount  and  peak  rate  of  runoff.  Given 
Data  only  daily  total  rainfall,  these  must  be  estimated  with  concep- 

tual methods  based  on  hydrologic  experience.  The  methods  used 
in  EPIC  (Williams  et  al . 1984)  for  estimation  based  on  daily 
rainfall  are  used  here  with  minor  modification. 
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The  EPIC  model  uses  the  Soil  Conservation  Service  (SCS)  Curve 
Number  (CN)  runoff -estimation  method,  with  some  important  modi- 
fications. The  basic  method  is  fully  explained  in  the  SCS 
Hydrology  Handbook  (USDA  SCS,  1972).  Briefly,  runoff  amount 
(Q)  is  related  conceptually  to  precipitation  amount  (P) , a 
soil-water  storage  value  (sy) , and  an  initial  abstraction  (Ia) 
as  follows: 


(P  - I )2 
a 

P + s - I 

[47] 

w a 

= 0.2s 

[48] 

w 

0 100  .. 

- 254(cn  • x> 

[49] 

The  method  gets  its  name  from  the  curve  number  parameter,  CN . 

As  indicated  in  equation  [49],  CN  may  vary  from  1 to  100.  P, 

Q,  I , and  s are  expressed  here  in  millimeters. 

^ a w r 

Equation  [47]  was  originally  described  (USDA  SCS,  1972)  as 
estimating  storm  runoff  amounts  rather  than  daily  amounts. 
Nevertheless,  it  is  applicable  just  as  well  (or  better)  for 
larger-period  rainfalls.  The  CN  is  conceived  as  a soil  parame- 
ter, with  a mean  value  for  average  soil  conditions.  A speci- 
fied method  modifies  the  curve  number  for  wet  or  for  dry 
initial  soil  conditions.  The  EPIC  modifications  used  in  the 
daily  option  of  Opus  relate  the  value  of  CN  to  current  soil- 
water  content  (Williams  et  al . 1984).  Opus  also  uses  the 
method  of  EPIC  to  estimate  peak  runoff  rates  for  runoff  events, 
based  on  estimated  time  of  concentration,  monthly  distribution 
of  rainfall  intensities,  and  ratio  of  maximum  30-minute  and 
24-hour  rainfall  depths  for  the  location  of  interest  (Williams 
et  al.  1984) . 

Unlike  EPIC,  Opus  provides  (and  the  authors  encourage  the  use 
of)  an  option  to  use  the  CN  rainfall/runoff  relation  described 
above  as  an  expected  value  rather  than  a prediction.  Thus  for 
long-term  simulation  sequences,  individual  runoff  amounts  are 
randomly  selected  from  a rational  distribution  about  this  mean. 
The  physical  limits  of  a runoff/rainfall  ratio  are  0 and  1,  and 
a two-parameter  beta  distribution  with  these  limits  is  used. 
Figure  18  illustrates  the  shape  of  the  curve  described  by 
equation  [47]  and  also  the  typical  scatter  of  runoff  data  from 
two  catchments  with  rather  long  records  (Montgomery  1980) . 

Some  scatter  of  CN-estimated  runoff  occurs  because  the  CN 
changes  with  initial  soil  water  content.  But  as  figure  19 
shows,  this  scatter  is  a minor  amount  compared  to  scatter  in 
recorded  data.  Scatter  about  such  a rainfall- runoff  curve  is 
to  be  expected  for  many  reasons  but  especially  because  of  the 
variation  in  storm  intensities.  This  variation  is  important 


Simulation  of  Runoff 
From  Breakpoint 
Rainfall  Data 


information  not  measured  by  daily  or  storm  rainfall  totals. 
Figure  20  illustrates  the  beta  function  distribution  of  runoff 
used  in  Opus  for  several  given  runoff/rainfall  ratios. 


Given  adequate  rainfall  data,  an  entirely  different  simulation 
method  is  available  with  Opus.  A deterministic,  distributed 
simulation  of  the  dynamic  response  of  soil  and  surface-water 
flows  is  done  when  breakpoint  or  pluviograph  data  are  provided. 


Infiltration 

Flow  within  the  soil  itself  is  assumed  to  be  described  by 
Darcy's  law  and  continuity  of  mass  as  expressed  in  equation 
[21].  Infiltration  is  flow  through  the  soil  surface,  and  the 
fundamental  relation  for  infiltration  is  also  based  on  Darcy's 
law  and  an  expression  of  continuity  for  flow  across  the  sur- 
face. Assuming  that  initial  water  content  9.  is  relatively 
uniform  in  the  surface  region,  and  with  infiltration  rate  f ( t ) , 
the  basic  relation  for  infiltration  (Smith  1981b)  is 

9 3 

1(f)  - f [50] 

**  9 . s 

i 

in  which 

D(0)  is  diffusivity,  defined  as  K(9)dip/d9  (mm2/min)  , 

I is  depth  of  infiltration  from  start  of  rainfall 

(mm) , 

f is  rate  of  infiltration  (mm/min) 

and  other  symbols  are  defined  in  chapter  4. 

Two  types  of  infiltration  conditions  may  be  encountered  at  the 
soil  surface.  The  first  is  an  imposed  ponding  condition,  such 
as  with  flood  irrigation.  This  condition  effectively  imposes  a 
fixed  soil -water  head  at  the  surface,  and  equation  [50] 
describes  infiltration  rate  (f)  for  all  values  of  I > 0 under 
these  conditions.  The  other  condition  is  the  common  rainfall, 
or  flux  boundary  condition.  For  this  case,  I at  first  increas- 
es as 


I 


t 


r( t)dt 
o 


[51] 


until  the  surface  becomes  saturated,  and  the  boundary  condition 
changes  to  a fixed  head  of  0.  Consequently,  at  that  time  and 
beyond,  the  infiltration  capacity  is  controlled  by  the  condi- 


41 


tions  near  the  soil  surface.  The  time  of  this  change  of  con- 
trol is  the  time  of  ponding  (t  ),  after  which  equation  [50] 
describes  the  decay  of  f with  ?ncreasing  I.  These  two  rela- 
tions are  illustrated  in  figure  21.  In  other  words,  at  t = t , 
equation  [50]  indicates  the  relation  between  rainfall  rate,  ^ 
initial  water  content  (0 .) , and  ponding  depth,  because  the 
lower  limit  of  the  validity  of  equation  [50]  is  reached  when  f 
- x and  I at  ponding  = I . An  analytic  function  is  obtained 
from  equation  [50]  by  miking  assumptions  on  the  functional 
form  of  D that  allow  integration  (Parlange  and  Smith  1976).  The 
function  used  in  Opus  is  that  developed  by  Smith  and  Parlange 
(1978),  which  assumes  exponential  behavior  of  K near  9 , so 
that  D(0)  is  approximately  proportional  to  dK/d0 , and  equation 
[50]  can  be  integrated  to  obtain 


1(f)  = G(0.)  In  [^-] 

s 


[52] 


Capillary  parameter  G(9.)  is  a coefficient  depending  on  initial 
conditions  and  is  related  to  the  integral  properties  of  the 
soil  capillary  characteristics.  For  most  initial  water  con- 
tents, it  can  be  represented  by 

G(9 . ) - H (9  - 9.)  [53] 

l c s i L J 


in  which  H is  a capillary  scale  parameter  and  is  relatively 
independent  of  initial  water  content  (Smith  1983a) . The  value 
of  H can  be  obtained  from  inf iltrometer  experiments,  or  from 
the  soil  characteristics  as 

1 r° 

Hc  “ K K^)dV>  [54] 

S ^ 
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Thus  a rough  estimate  of  H can  also  be  regressed  from  soil 
texture  as  was  done  for  unsaturated  soil  hydraulic  parameters 
(Rawls  et  al . 1983).  Equation  [52]  is  universal  in  the  sense 
that  it  applies  to  the  determination  of  ponding  occurrence  and 
to  the  rate  of  infiltration  after  ponding,  and  also  describes 
infiltration  from  initially  ponded  (flooded)  soil  conditions. 
Other  closely  related  "universal'1  relations  may  be  similarly 
obtained  by  making  different  assumptions  on  the  nature  of  D in 
the  integration  of  equation  [50]. 

When  inverted  as  an  f(I)  relation,  equation  [52]  becomes 


f v exp (I/G) 

s exp ( I/G) -1 


[55] 


A relation  for  f as  an  explicit  function  of  t that  is  accurate 
to  within  a very  few  percent  over  most  time  ranges  of  interest 
can  be  obtained  from  any  of  the  "universal"  relations  (Smith 
1981b).  In  particular,  for  equation  [55],  the  first  two  terms 
of  the  Taylor  expansion  of  exp(x)  may  be  used  to  solve  for  an 
expression  useful  in  simulation.  Equation  [55]  is  first  inte- 
grated with  f expressed  as  dl/dt,  to  obtain 

K t = I + G[exp( -I/G)  - 1]  [56] 

Then  using  the  Taylor  expansion  and  the  finite  differential 
expression  that  f -*•  AI/At  = [I(t+At)  - I(t)]/At,  we  obtain 

AI  = G JW  - 2KgAT/(beG)  - B [57] 

in  which 

B = (1-b  )/b  , and 

be  — exp(?/G)e 

This  expression  is  convenient  for  advancing  in  small  time  steps 
during  model  simulation,  with  maximum  errors  on  the  order  of  1 
percent.  Mean  time-step  inf iltration- rate  f is  taken  as  AI/At. 
Note  that  this  expression  is  asymptotically  correct  at  both 
small  and  large  times.  At  small  I , f is  proportional  to  1/1, 
or  I is  proportional  to  the  square  root  of  time.  As  I -+  «,  f 
approaches  Kg  asymptotically. 


Treating  Storm  and  Soil  Complications 

To  use  this  infiltration  theory  in  field  runoff  simulation,  one 
should  be  able  to  account  for  a variety  of  complexities  pre- 
sented in  nature.  This  section  describes  how  some  common 
complexities  are  treated  by  Opus.  It  should  be  kept  in  mind 
that  allowance  for  deterministic  complications  by  the  methods 
outlined  here  can  improve  simulation  accuracy  only  to  a certain 
degree.  The  major  ground  between  mathematical  simulation  and 
natural  system  performance  is  occupied  by  the  complexities  of 
spatial  and  temporal  variations  that  are  always  found  in  catch- 
ment characteristics.  Treatment  of  such  variability  is  within 
our  capability,  given  large  capabilities  of  computational  time 
and  storage,  but  the  collection  of  information  to  describe 
natural  spatial  variability  is  usually  quite  uneconomic.  In 
some  cases,  the  variability  of  nature  coupled  with  the  sensi- 
tivity and  nonlinearity  of  hydrologic  processes  can  cause  the 
variability  to  dominate  the  response  characteristics  of  a 
system.  This  possibility  should  not  be  used,  however,  to 
dismiss  the  value  of  physically  sound  simulation  models. 
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Rainfall  Hiatus.  One  common  complication  is  the  occurrence  of 
periods  within  a storm  in  which  rainfall  rate  falls  below  . 
Runoff  is  not  generated  during  these  periods,  and  any  runoff 
water  on  the  surface  from  previous  rainfall  will  remain  subject 
to  infiltration.  The  desaturation  of  the  surface  soil  proceeds 
in  accordance  with  Richards'  equation,  subject  to  complexities 
of  soil  structure  including  soil-water  hysteresis.  At  large 
times  after  the  beginning  of  such  a draining  period,  given  that 
r is  constant  and  less  than  K , the  surface-water  content  (8  ) 
approaches  the  kinematic  value  (Smith  1983b) , which  is 

*o  - "r  + V/Ks)-1/£(rs  - «r>  [581 

For  a rainfall  hiatus  in  Opus,  9 is  computed  from  r as  this 
equation  shows,  and  is  approached  asymptotically  with  a transi- 
tion half-life  of  approximately  30  minutes.  The  asymptotic 
approach  is  made  to  be  a function  of  the  proportion  of  surface 
still  covered  with  water,  the  amount  of  rainfall  already  infil- 
trated, and  the  hydraulic  conductivity  of  the  soil,  as  follows 
(Woolhiser  et  al . 1990): 

= flj-1  + [0  - 0j_1][l  - exp(-c  At)]  [39] 

in  which 

j is  an  index  for  successive  time  steps, 


At  is 


t . 
J 


cj-r 


c 

r 


is  a rate  parameter  estimated  as 


Cr  = KS/(I+0-1)  [60] 

This  is  a conceptual  relation  but  is  similar  to  one  that  has 
proved  successful  in  simulating  data  from  Walnut  Gulch  Experi- 
mental Watershed.  This  relation  provides  that  redistribution 
proceeds  more  slowly  for  large  rains  (large  I)  and  for  tighter 
soils  (small  K ) . 

This  approximation  of  surface-water  desaturation  is  applied 
only  during  the  course  of  a storm,  with  total  poststorm  redis- 
tribution calculated  by  Richards'  equation  as  described  above. 
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Infiltration  Through  a Surface  Crust  Layer.  Another  important 
complication  to  be  treated  in  natural  conditions  is  vertical 
anisotropy  in  the  soil  profile.  Opus  treats  a general  three- 
layer  case  in  the  simulation  of  infiltration  dynamics.  The 
surface  soil  is  assumed  to  lie  above  a subsoil  layer,  which  may 
be  assigned  properties  by  the  model  user  to  represent  a plowpan 
condition.  The  lower  soil  is  assumed  by  the  infiltration  cal- 
culations to  be  either  equal  to  or  less  than  the  surface  soil 
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in  hydraulic  conductivity.  This  assumption  is  not  required  for 
calculation  of  soil-water  movement  (equation  [26]).  In  addi- 
tion, agricultural  soils  are  assumed  to  form  surface  crust 
after  mechanical  disturbance,  due  to  the  impact  of  energy  of 
subsequent  rainfall.  The  following  method  of  treating  the 
effect  of  this  crust  is  not  theoretically  complete,  but  it  does 
provide  for  modification  of  infiltration  performance  in 
response  to  management  changes  of  the  surface -soil  conditions. 


The  crust  layer,  when  formed,  is  arbitrarily  assumed  to  com- 
prise the  top  10  mm  of  surface  soil.  Although  this  is  certain- 
ly not  accurate  for  all  cases  and  all  soils,  it  should  be  noted 
that  the  crust  thickness  and  the  crust  hydraulic  conductivity 
(expressible  jointly  as  a crust  resistance)  play  partly  inter- 
changeable roles  in  infiltration  dynamics  (Smith  1990) . 


The  Opus  infiltration  model  assumes  that  the  effective  saturat- 
ed hydraulic  conductivity  of  the  crust  (K  ) is  equal  to  or 
less  than  that  of  the  underlying  surface  soil  (K^)  . The  case 
of  Kg  > K ^ is  both  unusual  and  relatively  unimportant  in  its 
effect  on  infiltration.  When  the  ratio  of  r to  K is  large, 
ponding  could  occur  while  the  wetting  front  is  stitl  within  the 
crust  layer,  I < I . I is  I required  to  just  fill  the  crust 
layer.  In  this  case,  calculation  of  t is  straightforward, 
using  Kgc  in  equation  [56],  ^ 


Ponding  can  more  commonly  occur  after  the  wetting  front  has 
passed  the  crust  layer,  I > I . This  is  the  more  complicated 
case,  because  the  crust  and  underlying  soil  both  play  important 
roles  in  the  net  infiltration  dynamics  (Smith  1983b,  1990). 
Since  the  hydraulic  conductivity  of  the  crust  layer  has  been 
assumed  to  be  lower  than  that  of  the  underlying  soil,  the  final 
infiltration  rate,  or  effective  K , can  be  obtained  by  analysis 
of  steady  unsaturated  flow  across  isotropic  layers  (Smith 
1990) . 


The  vertical  distribution  of  matric  water  potential  in  the 
surface  region  will  generally  be  as  shown  in  figure  22.  At 
longer  times,  the  gradient  in  the  underlying  soil  becomes 
small,  the  value  of  becomes  closely  uniform,  and  the 
greater  gradient  of  ip  occurs  in  the  crust.  The  following 
equation  arises  from  matching  fluxes  into  and  out  of  the 
boundary  between  crust  and  subsoil,  assuming  that  head  at  the 
bottom  of  the  crust  (ip  ) equals  head  at  the  top  of  the  subsoil 
(V^)  (this  potential  is  termed  V^)  : 


1 


dtp 

dz 


MV>b) 


[61] 
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Since  the  right  side  of  equation  [61]  is  a constant  at  steady 
flow,  this  equation  applies  throughout  the  crust  layer  at  large 
times.  It  may  be  integrated  from  0 to  to  yield 


'^b d^> 

Kl(0b)/Kc(0)_1 


[62] 


The  function  K (V>)  is  described  in  equations  [2]  plus  [4],  and 

K1(^)  is  a single  value  satisfying  this  equation.  Equation 

[62]  is  integrated  by  trial  or  iteration  to  find  and  K1(t/>^), 

which  is  the  asymtotic  infiltration  rate  for  this  crusted 

profile  (K  ) . 
r se 

In  most  cases,  the  best  single  value  of  G to  use  when  the  wet- 
ting front  has  passed  into  the  subsoil  from  a crust  is  that  of 
the  subsoil  (Gx) . This  conclusion  comes  from  a study  of  infil- 
tration under  various  conditions  with  various  crust  properties 
using  solutions  of  equation  [26].  At  short  times,  equation 
[55]  is  accurately  approximated  by  the  expression 

f = KG/I  [63] 

This  indicates  that  soils  with  equal  products  K G(0.)  act  simi- 
larly at  short  times.  When  the  subsoil  has  a higher  K , and 
the  value  of  K G x is  lower  than  K G (using  K determined 
from  solution  of  equation  [62]),  thl  value  Gx  Ifiould  be  used  to 
simulate  infiltration  after  the  crust  has  been  passed.  When 
K Gt  is  larger  than  KG,  it  is  more  accurate  to  use  a weighted 
vafue  of  G for  I > I . S To  estimate  this  weighting,  G'  is  cal- 
culated, which  is  the  value  of  G such  that  G'K  = GK  at  I = 

I . Then  as  I increases  beyond  I , G for  use  m estimating 
ffl)  is  found  as  a weighted  average  of  G'  and  G:  (Smith  1990), 
as  follows: 

Ge  = (I-Ic)G'/I  + lQGl/I  [64] 

in  which  G is  the  effective  value  to  use  in  estimation. 
Although  this  procedure  may  seem  obscure,  it  is  handled  auto- 
matically in  Opus.  Validity  of  this  weighting  method  is  illus- 
trated for  one  case  in  figure  23. 


Crust  Development.  The  mechanics  of  transient  surface-crust 
layer  development  in  agricultural  conditions  is  not  understood 
sufficiently  to  allow  reliable  quantification  based  on  measur- 
able soil  parameters.  Crust  development  apparently  depends  on 
such  things  as  particle-size  distribution,  content  of  clay  and 
organic  matter,  exposure  to  rainfall  energy,  and  distribution 
of  rainfall  energy  (Chevalier  1984).  This  is  a subject  needing 
much  more  investigation.  It  is  known  that  the  hydraulic  con- 
ductivity of  crusts  developed  under  rainfall  is  commonly  an 


order  of  magnitude  less  than  that  of  the  parent  soil,  with  some 
noted  exceptions  (Chevalier  1984).  Further,  in  the  absence  of 
soil  cohesion,  the  crust  development  is  probably  limited  to 
simple  increases  in  bulk  density  due  to  packing. 

The  hypotheses  used  in  Opus  to  estimate  crust  formation  are 
illustrated  in  figure  24.  The  upper  graph  represents  a hypoth- 
esis for  the  effect  of  clay  content  on  the  ultimate  crust 
hydraulic  conductivity  that  can  be  formed  by  the  action  of 
rainfall.  The  lower  graph  illustrates  a hypothesis  for  the 
rate  of  crust  formation;  it  is  assumed  that  the  crust  formation 
is  an  exhaustion  function  related  to  the  accumulated  amount  of 
rainfall  energy- intensity  (El)  as  described  for  the  USLE  factor 
in  predicting  erosion.  The  storm  effective  El  used  in  crust- 
formation  estimates  is  directly  reduced  by  the  ratio  of  covered 
(plants  and  mulch)  surface  to  total  surface  of  the  catchment. 
Opus  uses  the  EPIC  regression  estimate  (Williams  et  al . 1984) 
to  obtain  an  estimate  of  El  when  daily  rainfall  data  are  used, 
and  finds  approximate  storm  El  using 

El  - 10  r2 

when  breakpoint  rainfall  rates  (r)  are  given  in  millimeters  per 
minute . 

Crusts  are  assumed  to  be  destroyed  by  cultivation  operations 
that  disturb  the  soil  to  any  depth.  Because  some  operations 
will  not  completely  destroy  an  existing  condition,  the  speci- 
fied efficiency  of  cultivation  (e  ) determines  the  starting 

state  (El  ),  as  follows; 
o 

El  - -4007^T  [65] 

In  part,  this  conceptually  represents  the  effect  that  cultiva- 
tion has  on  opening  the  soil,  which  goes  beyond  merely  destroy- 
ing a previously  formed  crust.  Thus  the  soil  hydraulic  conduc- 
tivity is  briefly  enhanced  (negative  El  ) as  a leftward  exten- 
sion of  the  lower  graph  in  figure  24.  Succeeding  rains  quickly 
increase  the  summed  El  into  the  positive  region  (e.g.,  clod 
dissolution) . 


Dynamics  of  Surface  Water 

When  the  infiltration  simulation  model  begins  to  produce 
rainfall  excess,  water  begins  to  move  over  the  surface  (except- 
ing the  case  of  surface  furrow  storage)  along  the  described 
topography  as  discussed  above.  Because  Opus  is  limited  to 
areas  containing  insignificant  spatial  variability  of  soil  con- 
ditions, all  parts  of  the  catchment  are  assumed  to  begin  to 
produce  runoff  at  the  same  time.  The  tracking  of  spatial  and 
temporal  variations  in  surface -water- flow  conditions  is  crucial 
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in  simulating  not  only  the  discharge  at  the  catchment  outlet, 
but  also  the  transport  processes  within  the  catchment  that  pro- 
duce the  pollutant  load  of  dissolved  and  suspended  material  at 
the  outlet. 

The  de  Saint-Venant  equations  for  surface  water  flow  have  prov- 
en to  be  a reasonably  accurate  and  efficient  means  to  describe 
the  dynamics  of  catchment  runoff.  Two  alternate  well-known 
approximations  are  used  in  Opus  to  solve  these  equations.  The 
kinematic  method  assumes  that  an  equation  of  mass  balance  and  a 
rating  equation  are  sufficient  for  all  but  very  flat  slopes  to 
describe  the  flow  of  surface  water.  The  rating  equation 
describes  a consistent  relation  between  discharge  and  flow 
depth  or  area  for  a given  section  geometry.  The  Manning  uni- 
form flow  relation  is  such  a rating  and  is  the  one  used  in  this 
model.  It  relates  flow  velocity  to  hydraulic  radius.  Hydrau- 
lic radius  is  defined  as  the  flow  cross-section  area  divided  by 
the  wetted  perimeter.  For  very  low  slopes,  a diffusive  wave 
equation  is  used,  slightly  more  complicated  than  the  kinematic 
wave  equation,  but  amenable  to  rapid  solution. 

Mass  balance  can  be  expressed  in  differential  form  as 


da  dq 
dt  dx 


r - f 


[66] 


in  which 

a is  cross-sectional  area  of  flow  (m2), 
x is  distance  from  beginning  of  runoff  path  (m) , 
q is  local  discharge  (m3/min) , 
t is  time  in  minutes, 

r is  current  rainfall  rate  (m/min) , and 
f is  local  current  infiltration  rate  (m/min) . 

The  Manning  equation  for  normal  flow  can  be  expressed  as 

q = 60  ^|PRb  [67] 

in  which 

S(x)  is  local  water  energy  slope, 
n is  Manning  roughness  coefficient, 

p is  wetted  perimeter  (m) , 

R is  hydraulic  radius  (m) , and 

b is  hydraulic  exponent,  = 5/3  for  Manning 

relation. 


Equations  [66]  and  [67]  are  also  tied  together  by  the  fact  that 
a = pR.  Thus,  equation  [66]  may  be  written  in  finite  differ- 
ence form  for  solution  as  follows: 


+ (1-w) [r-f ]| 


= F(R.)  = 0 


[68] 


in  which 

i refers  to  sequential  increments  along  the  x path, 
i=l  to  N, 

j refers  to  sequential  time  increments, 
w is  width  of  the  section  from  which  rainfall 
excess  contributes  (m) , 
u>  is  time  step  weight  (0  to  1)  , 

Ax  is  incremental  flow  length  = L/(N-1). 

When  the  flow  section  geometry  is  given,  the  relation  between 
a,  R,  and  p can  be  obtaine^  since  q^  = q(R^) . Thus  equation 
[68]  can  be  solved  for  R.~*  by  iteratively1making  F(R.)  ■=  0. 
The  Newton-Raphson  method  for  sets  of  equations  is  usei,  and 
the  matrix  equation  has  the  form 


J ( 3 , n)  5 (R)  = F(R) 


[69] 


where 


J(3,n)  is  the  tridiagonal  Jacobian  matrix,  containing 
terms  of  3F/3R, 

8 (R)  is  the  vector  of  estimated  corrections  to  R,  and 

F(R)  is  the  vector  of  residuals  as  in  equation  [68]. 

The  upstream  boundary  condition  and  the  downstream  boundary 
condition  depend  on  the  flow  type  as  discussed  below. 

Kinematic  Wave  Flow 

Morris  and  Woolhiser  (1980)  demonstrated  that  one  can  properly 
assume  kinematic  flow  for  larger  slopes  or  for  shallower 
depths.  Their  criterion  is  expressed  as  a dimensionless  number 
k , as  follows: 


o 


k = S L/h 


o o o 


[70] 
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in  which 


S is  bed  slope, 

L°  is  flow  surface  length  (m) , and 

h is  normal  flow  depth  (m) . 

o 

When  k is  greater  than  about  5 for  low  Froude  numbers,  the 
kinematic  wave  approximation  can  be  made.  That  is,  the  slope 
term  in  equation  [67]  can  be  assumed  to  be  the  bed  slope:  S(x) 
= S (x) . For  rainfall  on  a kinematic  slope  with  an  upstream 
divide,  the  upper  boundary  condition  is  simply  that  q(x=0)  or 
r(x=Q)  = 0.  For  cases  where  a small  channel  has  an  upstream 
input  (q  ) , the  upstream  boundary  can  be  specified  as  R(x=0)  = 
R(q  ) . °This  requires  an  algorithm  to  invert  the  q(R)  relation 
(equation  [67])  for  a given  flow  geometry.  The  kinematic  flow 
assumption  requires  no  downstream  boundary  condition.  The  for- 
mulation given  above  for  the  f inite - difference  solution  allows 
either  sequential  solution  of  equation  [68]  from  i=l  to  the 
downstream  boundary  (i=N)  or  a matrix  solution  for  all  i at 
once  as  in  equation  [69]. 


Diffusive  Wave  Flow 

When  k is  smaller  than  5,  the  kinematic -wave  approximation  is 
generally  inadequate,  and  the  energy  slope  may  be  expressed  as 


S = S (x) 
o 


dh 

dx 


[71] 


i 

! 

' 


Use  of  equation  [71]  in  equation  [68]  makes  the  formulation 
more  complicated  and  nonlinear,  so  that  sequential  solution  is 
not  possible  as  in  the  kinematic -wave  case.  In  Opus,  equation 
[71]  is  approximated  by  S - S (x)  - AR/Ax.  In  addition,  a 
relaxation  coefficient  is  required  in  iteration.  This  means 
that  stability  of  convergence  requires  that,  for  the  (k+l)th 
iteration , 

h+i  - «5IW  + h 

in  which  is  the  correction  to  R found  from  equation  [69], 

and  k is  the  iteration  number.  The  upstream  boundary  condi- 
tion for  runoff  from  a plane  beginning  at  a divide  or  for  a 
channel  with  no  upstream  inflow  is  S(0)  =0,  or  dh/dx(x=0)  = 
-So(0).  For  flow  with  upstream  inflow,  the  upper  boundary 
condition  at  x = 0 is  the  normal  flow  depth  for  that  inflow  and 
flow  geometry  (equation  [67]).  The  usual  downstream  boundary 
condition  diffusive  wave  flow  is  a critical  depth  condition, 
which  assumes  an  overfall  of  some  kind.  This  is  not  presumed 
to  be  universally  true,  but  it  is  not  a severe  assumption  and 
allows  convenient  computation  of  a flow-depth  relation.  The 
alternative  is  more  complex  backwater  calculations. 
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Estimation  of  Hydraulic  Roughness 

The  user  of  Opus  may  specify  a value  of  Manning  n to  represent 
roughness,  which  is  used  until  a surface-altering  management 
operation  takes  place.  For  natural  untilled  catchments  or  for 
simulating  other  fixed  surfaces,  this  value  applies  throughout 
the  simulation. 

When  the  field  is  tilled,  an  empirical  method  is  used  to  esti- 
mate the  resultant  value  of  n.  It  is  assumed  to  be  the  sum  of 
four  factors:  a bare  soil  value  for  a smooth,  clean  soil  sur- 
face; a roughness  caused  by  surface  residue;  a form  roughness 
caused  by  clods  and  surface  irregularities  formed  during  plow- 
ing; and  for  unfurrowed  surfaces,  the  roughness  effect  of  plant 
stems.  For  many  partially  cohesive  soils,  the  actual  form 
roughness  from  plowing  can  vary  considerably,  depending  on  soil 
wetness  at  plowing.  The  relationship  is  unknown  and  is  the 
subject  of  considerable  research.  The  relation  here  is  not 
presumed  to  be  universal. 


The  value  of  n for  bare  smooth  soil  is  assumed  to  be  0.02. 

Form  roughness  (f^)  varies  with  the  type  of  operation  and 
equipment.  A table  of  values  as  a guideline  is  given  in  the 
User  Manual  (vol.  II),  taken  partly  from  work  done  for  EPIC 
(Williams  et  al.  1984).  Values  range  from  6 to  50  and  repre- 
sent mean  height  of  microroughness  forms  produced  (in  mm).  The 
resulting  roughness  effect,  in  terms  of  a partial  Manning  n 
(n^.)  , is  estimated  (G.R.  Foster,  personal  communication)  as: 

nf  =■  0 . 00087ff1 ' 25  [72] 


The  roughness  resulting  from  surface  residue  is  estimated  fol- 
lowing the  work  of  Foster  et  al . (1982).  It  is  a function  of 
effective  residue  weight  (m  ) (kg/ha)  and  a coefficient 
(c  ) that  changes  betweenrfurrowed  and  unfurrowed  flows.  The 

efflctive  mulch  cover  (F  ) is  first  estimated  as 

m 

F = 1 - exp ( -m  /4428)  [73] 

m r res'  L J 


The  partial  residue  roughness  (n  ) is  then  estimated  accord- 
ing to  the  data  from  Foster  et  aI?S(1982)  by 


n 


res 


= 1/[1  + 


( c /m  ) 2 
res'  res 


0.61 


[74] 


The  value  of  c is  10.67  for  furrowed  surfaces  and  3.43  for 

1T6  S 

unfurrowed  surfaces.  Finally,  the  weighted  sum  of  form  and 
mulch  factors  becomes  a partial  n value  represented  by 


n 


rf 


F n + 
m res 


^1Fm^nf 


[75] 
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The  net  n also  includes  a factor  for  estimated  stem  density 
when  the  surface  is  unfurrowed  and  flow  must  pass  by  plant 
stems.  This  is  a function  of  the  plant  mass  per  unit  area  of 
surface.  The  final  estimating  equation  becomes 


n = 


n 


rf 


1.5  1.5  0.667  „ , , . 

+ *b  ] + 3-6E-6Miv/yT 


[76] 


in  which  M 
height  of 


lv  1S 
plant 


mass  of  plant  above  ground  (in  kg/ha)  and  y 
(in  m) . 


P 


is 


Flow  Geometry  Relations 

The  flow  section  geometry  for  which  equation  [68]  can  be 
applied  is  quite  arbitrary,  including  plane  flow;  furrows;  and 
trapezoidal,  triangular,  or  rectangular  channels.  When  a fur- 
row or  channel  is  treated,  often  a relatively  flat  bottom  is 
formed  by  local  transport  and/or  deposition  of  sediment.  This 
surface  can  have  a very  different  roughness  than  do  the  sides 
of  the  section.  In  such  a case,  the  local  net  value  of  Manning 
roughness  (n)  may  be  formed  as  a weighted  sum  of  the  roughness 
(n)  of  each  part  of  the  wetted  perimeter,  as  follows: 

RV  + "zV’M 0,5 

in  which  subscript  b refers  to  the  bottom  section,  and  sub- 
script z refers  to  side  slopes  of  the  section. 


Time  Steps 

For  surface-water  routing  with  this  scheme,  time  steps  are  the 
same  as  those  for  the  infiltration  calculations,  described 
above.  Hydraulic  conditions  on  the  catchment  are  checked  at 
each  time  step,  and  time  intervals  are  reduced  if  velocities 
between  x increments  are  large  enough  to  cause  potential  insta- 
bilities. These  reduced  time  steps  are  then  used  in  the  next 
time  increment  for  infiltration  simulation.  The  largest  time 
steps  allowed  are  those  read  from  the  breakpoint  rainfall 
record. 


Irrigation  Flows 

The  numerical  scheme  described  above  is  applicable  for  estimat- 
ing (a)  the  advance  of  water  down  a furrow  from  border  irriga- 
tion or  (b)  any  runoff  produced  by  sprinkler  irrigation. 
Sprinkler  irrigations  are  treated  like  rainfall  inputs,  except 
that  they  are  assumed  to  be  controlled  so  that  runoff  is  not 
produced.  For  border  irrigation,  the  limits  of  practicality  do 
not  allow  Opus  to  divide  long  furrows  into  enough  distance 
increments  to  produce  a precise  prediction  of  flow  advance,  but 
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Erosion  and 
Transport  of 
Sediment 


Opus  provides  a reasonable  numerical  approximation.  The  prin- 
cipal additional  feature  required  for  irrigation- advance 
simulation  is  the  proper  calculation  of  the  spatial  variation 
of  infiltration  rate,  which  depends  on  the  time  of  initial 
wetting  of  a given  point  along  the  furrow.  Thus,  infiltration 
rate  (f)  and  depth  (I)  are  separately  calculated  for  each 
increment  node  (x.)  along  the  plane.  Net  infiltration  for 
water-balance  calculations  is  obtained  by  summing  over  all 
locations  along  the  field. 


The  methods  used  by  Opus  to  estimate  sediment  production 
resulting  from  rain  on  a catchment  correspond  in  approach  to 
the  methods  available  for  estimating  runoff.  The  sediment 
option  accompanying  the  daily  runoff  method  is  a lumped- 
parameter  expression  that  corresponds  in  sophistication  to  the 
SCS  Curve  Number  runoff  method.  The  corresponding  option 
associated  with  breakpoint  rainfall  is  a spatially  distributed 
treatment  of  sediment  transport  similar  to  KINEROS  (Smith 
1981a) . The  option  has  been  expanded  to  allow  simulation  of 
sediment  with  up  to  five  representative  classes  of  particle 
size . 

For  all  cases,  information  on  the  rate  of  erosion  (or  deposi- 
tion) as  a function  of  soil,  surface,  cover,  and  local  hydrau- 
lic conditions  is  crucial  for  accurate  estimation  of  sediment 
transport.  Opus  makes  use  of  information  on  the  local  condi- 
tions that  is  available  from  other  parts  of  the  model;  included 
is  information  on  plant  and  mulch  cover,  and  soil -surface 
conditions . 


Sediment  Production  for  Daily  Simulation 


The  daily  amount  of  sediment  that  is  produced  by  the  daily  run- 
off is  estimated  with  a variation  of  the  well-known  USLE 
(Wischmeier  and  Smith  1978) , which  is  MUSLE  (Williams  1975) . 
USLE  is  a lumped-parameter  method  of  estimating  net  erosion  by 
use  of  a log- regression  expression  involving  a parameter  esti- 
mate of  each  of  several  major  erosion  factors,  as  follows: 


Q = R K L S 4>  P 


u u u u 


u 


in  which 


[78] 


Q 

R 

K 

L 

S 

4> 


s 

u 

u 

u 

u 


p 

u 


is  net  storm  or  daily  soil  loss  (kg/m2), 

is  storm  or  daily  erosivity  (Newtons/hr) , 

is  soil  erodibility  factor  (kg'hr/Newton/m2) , 

is  slope  length  factor, 

is  slope  steepness  factor, 

is  coefficient  for  cover  and  management, 

is  a factor  for  effects  of  supporting  management 

practices . 
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Use  of  USLE  and  detailed  discussion  of  the  coefficients  and 
their  values  is  provided  in  a USDA  Handbook  on  the  subject 
(Wischmeier  and  Smith  1978).  In  MUSLE,  R is  modified  from  the 
original  USLE  to  account  for  both  rainfall  erosivity  and  runoff 
erosivity.  The  estimating  equation  is  (Williams  and  Berndt 
1977) : 

R = 90 . 5 (Q ' q )°'56  [79] 

u np 

in  which 

Q is  runoff  volume  (m3),  and 
q^  is  peak  runoff  rate  (m3/s) . 

Clearly  this  is  a lumped  approximation,  but  the  method  is  test- 
ed and  has  complexity  appropriate  to  match  the  daily  hydrology 
estimation.  The  effects  of  impoundments  such  as  tile  outlet 
terraces  or  holding  ponds  are  simulated  by  modifications  to  the 
parameter  (P  ) . 


Distributed  Erosion  and  Sediment  Transport 

When  breakpoint  rainfall  is  used,  spatial  and  temporal  distri- 
bution of  runoff  hydraulics  is  simulated  and  can  be  used  to 
estimate  spatial  and  temporal  erosion  and  sediment  transport. 
The  general  governing  equation  for  sediment  transport  (Bennett 
1974)  in  time  (t)  and  along  distance  (x)  can  be  written  as 

k <aCsk)  + k <qCsk>  - d<x-C)k  - qs(x,t)k  1801 

in  which 

C is  sediment  concentration, 

S 

a is  cross-sectional  area  of  flow  (m2 ) , 

q is  water  discharge  per  unit  width  (m2/min) , 

d is  rate  of  erosion  or  deposition  at  the  bed  (m/min) , 

qg  is  local  input  of  sediment  (m3/m2/min) , and 

k"  is  a subscript  referring  to  a particle  size  class. 

This  equation  may  be  applied  and  solved  through  time  and  space 
in  conjunction  with  the  information  on  hydraulic  conditions  in 
time  and  space.  It  is  solved  in  conjunction  with  the  numerical 
solution  for  runoff  water.  If  the  user  knows  the  hydraulic 
variables  for  a time  step,  the  change  in  concentration  with 
time  over  space  can  be  found  directly. 


Surface  Erosion  and  Sediment  Transport 

Whether  the  surface  is  furrowed  or  unfurrowed,  flow  can  be 
thought  of  as  composed  of  both  "sheet"  and  rill  flow.  Erosion 
is  caused  by  the  shear  of  flowing  water  and  by  the  disturbance 
of  surface  soil  by  raindrop  impact.  For  this  case,  Foster 
(1982),  based  on  existing  experimental  data,  proposed  an  equa- 
tion for  the  detachment  of  soil  particles  by  rainfall  as 
follows : 

d.  = 0.0138  K r 2 [ 2 . 96 (S . )° ' 79  + 0.56]  <f> . [81] 

l u e 1 l J 1 J 

where 

r is  effective  rain  intensity  (mm/hr) , 

S.  is  sine  of  the  interrill  surface  - slope  angle,  and 

<f> ^ is  the  interrill  soil/crop  coefficient  as  follows: 

*i  = <t>u  (l-Fs)exp[-0. 21(hc/hb  - 1) 1 ‘ i0]  [82] 

where 

F is  fractional  total  cover  (plant  and  mulch)  over  soil, 

h is  water  flow  depth  for  given  cover  conditions  (m) , 

h^  is  water  flow  depth  for  bare  soil  (m) , 

<b  is  a surface-soil  residue  factor  defined  as 

u 

4>  = B (c  + c,t  )exp(-120f  ) [83] 

^u  s a b m r rs  L J 

in  which 


B is  a soil  consolidation  factor, 

s 

c and  c,  are  coefficients, 

cl  U 

t is  time  since  last  tillage,  and 

f is  residue  concentration  in  surface  soil  (10-100  mm) 
rS  (kg/m2/mm) . 

Foster's  concepts  grew  out  of  the  USLE  to  some  extent,  and  the 
<t>  factors  are  related  in  concept  to  the  <f>  of  the  USLE,  insofar 
as  they  represent  effects  of  management  on  the  erosion  process 
The  coefficients  are  discussed  in  Foster  et  al.  (1983). 

Interactions  of  rill  erosion  and  transport  capacity  are  calcu- 
lated with  a linear  deficit- reaction  relationship  as  follows: 

drk  - <dpr/gc,[gc  • gsL  1841 


55 


in  which 


is  actual  local  rill -detachment  rate  (kg/m2/min) , 
is  potential  rill  detachment  rate  (kg/m2/min) , 
is  transport  capacity  for  particle  class  k 
(kg/m2/min) 

is  current  transport  rate  for  particle  class  k, 
(kg/m2 /min) , and 

is  an  index  for  each  particle  - size  class. 


For  this  option,  Foster  et  al . (1983)  proposed  d to  be 
related  to  local  hydraulic  shear,  as  follows: 


d 

Pr 


[85] 


and  he  also  proposed  that  transport  capacity  is  similarly 
related  to  shear,  as  follows: 


g 


c 


3/2 

r 

s 


[86] 


in  which  parameters  bn  and  b are  defined  in  the  following 
(Foster  et  al . 1983): 


= 139K 

<t>  <t> 

[87] 

u 

u r 

= 188  - 

468f  + 907f  2 

; f 

< 0.22 

[88a] 

c c 

c 

= 130 

'fc 

> 0.22 

[88b] 

in  which 


f is  clay  fraction  in  soil,  and 

4>l  = exp[-1.8Mm]  [89] 

which  is  similar  to  the  factor  <f>  . M is  mulch-cover  density 
/ • i / 2 \ u m J 

(m  kg/m2)  . 


With  the  shear  terms  in  equations  [85]  and  [86]  dividing  out  in 
equation  [84]  , d /g  becomes  a constant  related  to  factors  <f> 
and  clay  content . ^rThus  equation  [84]  is  very  similar  to  the 
equation  used  in  KINEROS  (Woolhiser  et  al . 1990)  except  that 
variations  in  the  coefficient  are  still  found  from  changes  in 
the  <f>  factors  as  shown  in  equations  [87]  and  [88]. 

Potential  transport  capacity  is  found  using  the  transport 
capacity  relation  of  Engelund  and  Hansen  (1967)  as  follows: 

0.05  S u u^3 

c = ~ o- u 7 FT?  [90] 

smx  g2  h yd  (Ps*l)2  1 J 
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in  which 


u = flow  velocity  (ra/min) , 
u^,  = shear  velocity  (m/min)  , 

y^  =■  mean  particle  size  (m)  , 

= particle  specific  gravity, 
g = acceleration  of  gravity  (m/min2), 

h = flow  depth  (m) , and 

Csmx  = transPort  capacity  concentration  for  given 
particle  and  hydraulic  conditions. 

This  relation  is  used  because  it  has  shown  itself  to  be  as 
robust  as  any  for  shallow  flows  and  high  slopes,  which  may 
often  be  encountered  (Alonso  et  al . 1981).  It  has  been  dis- 
cussed thoroughly  elsewhere  (Alonso  1978). 


Concentrated  Flow  Erosion  and  Sediment  Transport 

Even  in  catchments  of  relatively  simple  geometry,  distributed 
surface  flows  soon  become  concentrated  into  rivulets  or  small 
channels.  On  managed  areas  these  channels  may  be  the  furrows 
themselves,  or  they  may  be  formed  when  flow  in  furrows  finds  a 
local  swale  and  forms  a cross-furrow  flow.  These  cases  often 
result  in  ephemeral  gullies  in  areas  where  erosion  potential  is 
high.  They  are  ephemeral  in  cultivated  lands  because  they  are 
obliterated  by  mechanical  cultivation.  Conversely,  these 
swales  may  be  places  for  upland  eroded  material  to  deposit. 
Foster  et  al . (1983)  suggested  that  erosion  rates  in  small 
channels  are  best  modeled  by  a relation  based  on  the  hydraulic 
bed  shear  and  including  consideration  of  a critical  or  thresh- 
old value,  as  follows: 

D,.  = c-.(1.35r  - t ) [91] 

f f s cr 

in  which 

Df  is  local  concentrated  flow- detachment  rate  (in 
kg/m2 /min) , 

c^  is  a coefficient, 

r is  hydraulic  shear  on  the  wetted  perimeter  of 
flow  (Newtons/m2),  and 

r is  critical  shear  above  which  erosion  occurs, 

C IT 

(in  Newtons/m2). 

This  value  is  considered  a potential  detachment  rate  and  is 
used  in  equation  [84]  along  with  transport  capacity  for  dis- 
tributed flow  to  determine  actual  local  erosion  (or  deposition) 
rate.  The  value  c^  is  suggested  by  Foster  et  al . (1980b)  to  be 
estimable  as  proportional  to  USLE  erodibility  (k^) ; in  Opus,  c^ 

is  taken  to  be  0.246*k  (metric  units). 

u 
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Erosion  and  Deposition  With  Mixed  Particle  Sizes 

The  transport  of  sediment  particles  is  possible  because  of  the 
dissipation  of  hydraulic  energy  from  the  flowing  water.  In  the 
erosion  and  transport  of  soils  and  sediments  of  a mixture  or 
range  of  particle  sizes,  different  particle  sizes  will  exhibit 
very  different  transport  capacities  (see  equation  [90]). 

Although  there  are  yet  no  definitive,  experimentally  based  laws 
governing  the  apportionment  of  transport  energy  among  particle  1 
sizes  of  a mixture,  some  reasonable  assumptions  can  make  possi-  j 
ble  a useful  model  for  the  process.  Opus  treats  a distribution 
of  particle  sizes  as  did  CREAMS,  that  is,  by  using  up  to  five 
different  classes  of  particle  size.  The  methods  of  Foster  et 
al.  (1985),  for  estimating  the  distribution  of  classes  in  the 
absence  of  specific  information,  are  available  as  a default  in 
Opus . 

. 

In  the  distribution  of  transport  capacity  among  sediment-size 
mixtures,  three  general  cases  exist.  In  the  first  case,  the 
transport  capacity  of  flow  may  exceed  the  current  sediment  load 
of  flow  for  all  particle  sizes  involved.  Here  it  is  reasonable 
to  assume  that  erosion  will  act  on  all  the  particles  of  the 
soil,  so  that  particle  sizes  are  eroded  in  direct  proportion  to 
their  availability. 

In  the  second  possible  case,  the  transport  capacity  of  flow  is 
inadequate  to  transport  the  sediment  load  for  any  particle 
sizes  involved.  For  this  case,  deposition  of  all  sizes  will 
occur . 

In  the  third  case,  conceptual  assumptions  are  required,  and 
present  erosion  theory  is  incomplete.  Here  the  transport 
capacity  may  be  greater  than  the  sediment  load  if  composed  of 
only  the  smallest  particle  class,  but  the  reverse  is  true  for  a 
load  composed  entirely  of  the  largest  particle  size.  The  model 
must  distribute  the  hydraulic  energy  among  the  particle  sizes 
that  make  up  the  actual  sediment  load.  In  Opus,  the  total 
transport  capacity  is  assumed  to  be  a weighted  sum  of  the 
transport  capacities  calculated  for  each  class  of  particle 
size.  The  particle-size-class  capacity  weightings  are  assumed 
to  be  inversely  proportional  to  the  fall  velocities  of  the 
class  of  particle  size.  Thus,  a particle  size  with  relatively 
high  transportability  will  tend  to  be  enriched  in  the  total  mix 
of  particles.  Since  smaller  particles  almost  always  have  a 
smaller  fall  velocity  relative  to  their  proportion  of  the  soil, 
they  tend  to  be  always  enriched  in  the  sediment  in  runoff. 

This  occurs  not  by  excessive  erosion,  but  by  preferential  depo- 
sition of  heavier  particles. 

In  determining  the  local  conditions  for  transport  of  several 
particle-size  classes,  the  total  transport  capacity  is  compared 
with  the  total  sediment  load.  If  excess  capacity  exists,  a net 
erosion  is  indicated,  but  Opus  assumes  that  erosion  must  take 


all  particles  in  proportion  to  their  distribution  in  the  soil, 
or  else  erosion  would  be  self - limiting  by  surface  armoring. 

Thus  erosion  is  distributed  according  to  the  proportion  of  par- 
ticle sizes,  with  a rate  necessarily  controlled  by  the  rate 
determined  for  the  largest  (and  least  transportable)  size.  On 
the  other  hand,  when  there  is  an  excess  of  total  load  over 
total  transport  capacity,  selective  deposition  may  take  place. 
In  this  case,  deposition  will  occur  for  the  larger  particle 
sizes,  to  the  extent  that  total  load  exceeds  total  capacity. 

It  is  possible  for  net  hydraulic  erosion  to  be  zero  but  with 
excess  of  capacity  for  small  sizes  and  deficit  of  capacity  for 
large  sizes.  This  is  conceptually  analogous  to  armoring  for 
riverbeds,  but  the  rapid  changes  in  shallow  flow  from  rainfall 
runoff  make  actual  armoring  unlikely  as  an  extensive  or  signif- 
icant watershed  phenomenon.  However,  selective  deposition  and 
consequent  enrichment  of  fines  is  common  and  is  simulated  by 
Opus  as  described  above . 


Evolution  of  Flow  Sections  With  Erosion  and  Deposition 

Erosion  and  deposition  often  occur  alternately  during  a runoff 
event  at  a given  location,  and  will  alter  the  geometry  of  a 
furrow  formed  by  a farming  implement.  Also,  small  channels  in 
a catchment  will  tend  to  have  a geometry  that  is  characteristic 
of  the  balance  of  erosive  and  depositional  forces.  With  some 
simplifying  assumptions,  one  can  estimate  the  direction  that 
will  be  taken  in  erosive  or  depositional  changes.  Opus 
includes  an  algorithm  to  estimate  the  changes  in  a furrow  or 
channel  section  that  result  from  the  loss  or  gain  of  a certain 
eroded  or  deposited  volume  (expressed  as  m3/m) . 

It  is  assumed  that  the  channel  cross  section  starts  as  a simple 
trapezoid  formed  as  a furrow  or  as  the  intersection  of  two 
slopes,  and  is  reset  by  mechanical  operations.  The  width  of 
the  bed  is  assumed  to  adjust  itself  until  critical  shear  is 
exceeded  only  on  the  bottom  of  the  section.  Changes  of  section 
shape  are  not  computed  during  a runoff  event,  but  are  effected 
after  the  event  when  flow  means  and  erosion  or  deposition 
totals  are  known.  Erosion  of  a trapezoidal  section  can  be 
accompanied  by  either  widening  or  narrowing  of  the  bottom 
section.  The  same  is  true  with  deposition  or  reduction  in 
section  area. 

In  CREAMS,  Foster  et  al . (1980a)  assumed  a simple  distribution 
of  the  local  shear  force  along  the  wetted  perimeter  and  devel- 
oped an  algebraic  method  to  predict  bottom  width  based  on  that 
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assumption.  If  the  flow  section  is  assumed  to  be  rectangular, 
the  hydraulic  radius  is  related  to  the  section  factor  as 
follows : 

R - [SF]3/8c  [92] 

m 

in  which 

SF  is  a section  factor,  Qn/7S  , and 

c^  is  a coefficient  that  depends  on  the  section  geometry. 

The  local  shear  stress  is  assumed  to  be  distributed  along  the 
wetted  perimeter,  increasing  from  zero  at  the  water  edge  to  a 
maximum  at  the  middle  of  the  bottom.  To  describe  this  distri- 
bution we  define  a normalized  distance  along  the  wetted  perime- 
ter (here  called  y^) , which  is  obtained  by  dividing  the  dis- 
tance along  the  perimeter  by  the  wetted  perimeter  (p) . The 
local  shear  stress  can  be  expressed  in  terms  of  critical  shear 
stress  as 

r = r B 
c 

with 


r = critical  shear  stress  - 7S  R 
c o 

in  which 


7 is  specific  weight  of  water,  and 
B is  a weighting  function  found  so  that 
rP/2 


1.0  = 


Bdx 


J0 


[93] 


These  equations  merely  express  the  concept  that  shear  stress  is 
distributed  along  the  distance  from  the  water  edge  to  the  mid- 
dle of  the  rectangular  channel  (y  = 0 to  p/2,  or  yA  = 0 to  y^  = 
0.5),  and  weighted  in  such  a way  that  the  total  shear  stress  is 
still  expressed  by  equation  [93].  If  c and  R in  equation  [92] 
are  expressed  in  terms  of  y and  combined  with  equation  [93],  a 

combined  factor  b can  be  obtained  as  follows: 
c 

bc(y*)  = [ SF] 3/8  <7S Q/rc)  [94] 

The  terms  on  the  right  can  be  obtained  independently  of  the 
section  proportions,  and  the  function  b (y  ) is  inverted  so 
that  y^  is  found  as  the  dimensionless  depth  of  the  rectangle 
with  r at  the  edge  of  the  bottom  just  equal  to  r . The  bottom 
width  immediately  follows  from  the  rectangular- shape  assumption 
and  can  be  expressed  as 

wf  = [ SF(  1 - 2h)  ] 3/8  y*'V8  [95] 
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This  provides  a means  to  estimate  the  width  of  the  bottom  when 
the  shape  of  the  section  is  being  altered  by  scour  or  deposi- 
tion. This  method  is  used  in  Opus  with  sloping  sides  rather 
than  rectangular  sections.  Figure  25  illustrates  how  the 
section  is  assumed  to  change.  At  the  end  of  a runoff  event, 
the  amount  eroded  or  deposited  at  each  section  has  been  accumu- 
lated and  expressed  in  terms  of  a cross-sectional  area.  This 
is  used  along  with  the  width  found  from  equations  [94]  and  [95] 
to  solve  for  the  required  cross  section.  Given  a change  in 
section  area  and  new  bottom  width,  the  change  in  bottom  eleva- 
tion is  the  single  unknown  dimension. 


Impoundment  Storage  One  significant  management  practice  that  may  often  be  used  in 

controlling  pollution  from  agricultural  catchments  is  the 
inclusion  of  an  impoundment  element  in  the  field.  This  may 
take  the  form  of  a tile  outlet  terrace  or  a series  of  them,  or 
it  may  be  as  simple  as  a pond  at  the  lower  end  of  the  field. 

Of  course,  ponds  may  also  be  natural  parts  of  unmanaged  catch- 
ments. For  experimental  catchments,  a pond  may  be  formed  as 
backwater  to  a measuring  weir  or  flume. 

The  option  for  daily  rainfall  sediment  (described  above)  simu- 
lates the  effects  of  a pond  on  sediment  and  sediment-associated 
pollution  by  modification  of  the  USLE  "P"  factor  (P  ) , as 
specified  in  USLE  methodology.  The  USLE  handbook  describes  the 
appropriate  P^  factor  changes  for  pond  effects. 


Routing  of  Water  Through  a Pond 

Given  the  time  distribution  of  inflows  of  water  and  sediment  to 
a pond,  the  hydrologic  effects  of  the  pond  on  both  hydrograph 
and  sediment- concentration  distribution  can  be  analytically 
simulated.  This  is  done  in  the  option  for  breakpoint  rainfall 
hydrology  and  sediment. 

The  pond  is  described  for  purposes  of  simulation  either  by  an 
equation  for  depth-area  relationship,  or  by  the  slope  geometry 
of  the  surfaces  converging  at  the  pond  location.  That  is,  the 
description  of  the  slope  of  two  planes  and  the  slope  of  their 
intersection  (presumably  a channel  to  which  they  contribute, 
plus  the  slope  of  the  face  of  a dam  across  the  channel)  can  be 
converted  into  an  equation  for  the  relation  between  depth  of 
water  and  surface  area  contained  in  this  inverted-pyramid- 
shaped  "basin." 

Either  method  produces  an  equation  with  three  parameters  relat- 
ing pond  water  surface  area  (A)  to  depth  (h)  as  follows: 

c 

A(h)  = A + b h P [96] 

P P 
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in  which  A 
pond. 


and  0^  are  the  descriptor  parameters  for  the 


The  other  information  needed  to  route  water  and  sediment 
through  the  pond  is  the  outlet-rating  relation.  This  is  an 
equation  relating  discharge  from  the  pond  to  depth,  and 
includes  a threshold  depth  below  which  the  discharge  is  zero. 
The  pond  is  assumed  to  seep  water  into  the  soil  at  all  depths, 
using  the  saturated  hydraulic  conductivity  of  the  surface-soil 
layer  or  other  user- specif ied  limiting  conductivity  as  the  rate 
of  seepage. 

The  outflow  rating  may  be  given  either  as  an  explicit  relation 
of  discharge  to  depth  or  in  the  form  of  an  orifice  coefficient, 
which  is  translated  by  the  program  into  an  explicit  relation. 
The  general  rating  relation  for  outlet  discharge  is 


Q = c (h-h  ) 
o o z 


[97; 


in  which 


Q is  outlet  discharge  (m3/min) , 

h is  pond  depth  (m) , 

h is  the  threshold  depth  for  outflow  (m) , and 

7. 

cq,  d^  are  parameters. 

The  routing  of  runoff  through  the  pond  is  simply  the  solution 
of  the  differential  equation  for  a linear  storage  system,  as 
follows : 


dV 

dt 


m 


in  which 


[98] 


Q.  is  pond  inflow  discharge  (m3/min) , 
V is  storage  in  pond  (m3). 


Q^n(t)  is  known,  and  the  solution  is  obtained  through  the  rela- 
tion between  Q (t)  and  h of  equation  [97],  plus  the  relation 
between  the  pond  volume  V and  depth  y(t),  which  comes  from  the 
following  integration  of  the  area-depth  relation  equation  [96]: 


V(h) 


A(h)dh  = A h + 

0 p 


b 

P 

c +1 
P 


(c+1) 
h p 


[99] 
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In  Opus,  the  solution  is  obtained  on  successive  small  time 
steps  (At)  during  which  Q.  is  assumed  constant.  Using  dV(h)  = 
A(h)dh,  the  solution  is  actually  obtained  in  the  form 


dh  Q.  - Q (h) 
= m o 

dt  A(h) 


[100] 


Pond  storage  (except  for  the  smallest  ponds)  will  often  have  a 
considerable  damping  effect  on  the  peak  discharge  from  a small 
catchment.  A significant  value  of  h indicates  that  runoff 
from  small  flows  can  be  trapped,  and  pondages  in  general  can 
trap  most  sediment  carried  by  all  but  the  largest  flows. 


Routing  of  Sediment  Through  a Pond 


Each  class  of  particle  size  from  the  sediment- transport  simula- 
tion is  considered  separately  in  routing  sediment  through  a 
pond.  A given  sediment  particle  entering  the  pond  is  assumed 
to  settle  at  its  fall  velocity.  The  general  equation  for  mass 
balance  of  the  sediment  particle  class  j,  with  concentration  C. 
and  settling  velocity  v . , can  be  written  as  ^ 


d(C.V) 

— nr1 — = C..Q.  - C.Q  (h)  - A(h)C . v , 

dt  J ixin  jHov  J sj 


[101] 


We  assume  (as  above  for  pond  depth)  solution  during  a 
time  step  in  which  Q and  Q.  are  locally  constant, 
that  dV/dt  is  Q.  - Q , the  general  equation  in  terms 
becomes 


small 

Noting 

of  C. 

J 


C..Q. 

1 i in 


C . (Q.  + Av  ) 

1 ^in  s 


V(t) 


[102] 


in  which  C..  is  input  value  of  C..  Equation  [102]  is  solved  by 
a transformation  of  variables,  aAd  the  solution  takes  one  of 
three  forms,  depending  on  the  relative  size  of  various  variable 
terms . The  forms  are 


Form  1 : For  relatively  constant  V,  or  Qin  = Qq, 

solution  for  C.  after  an  interval  At  is 
j 

C.  = C.  + (C.  - C.  )exp(-bAt/V) 

J jx  jo  jx' 


the 


[103] 


in  which 


C.  is  C.  at  time  0, 
J°  J 
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b 

is  Q.  + Av  . , and 
m sj 

I 

C. 

jx 

is  C . . Q . /b 

J i m' 

Form 

2 : For  |b/AQ|  greater  than  10, 

C.  = 
J 

0hx  + (CjQ  - Cjx)V/(bAt)  [ l-exp( -bAt/V)  ] 

[104] 

Form  3 : For  |b/AQ|  less  than  10, 

C.  = C.  + (C.  - C.  )[V  /(V  +AQAt)]b/AQ  [105] 

J jx  jo  jx  o'  o 

The  value  AQ  is  Q.  - Q . Time  steps  At  are  those  chosen  as 
appropriate  for  tiie  solution  of  equation  [100]  above. 

From  the  point  of  view  of  the  various  flow  paths  across  the 
impoundment,  given  a variety  of  flow-entrance  locations  and 
outlet  location,  the  pond  is  actually  a spatially  lumped  ele- 
ment in  the  method  as  described. 

After  the  cessation  of  input  to  the  pond,  its  settling  and  out- 
flow are  simulated  until  outflow  ceases,  with  time  steps  chosen 
to  optimize  the  computations.  The  fact  that  the  pond  is  on 
soil  with  a finite  loss  rate  assures  that  outflow  returns  to 
zero  at  a finite  time. 


Simulation  of  Snow  Precipitation  records  very  rarely  include  data  on  whether  the 
Budget  precipitation  in  winter  periods  occurs  in  the  form  of  snow  or 

rain.  Further,  since  Opus  is  designed  to  not  require  detailed 
daily  records  of  temperature,  there  is  no  way  to  accurately 
predict  the  occurrence  of  snow.  Even  when  records  of  daily 
maximum  and  minimum  temperatures  accompany  the  precipitation 
record,  often  it  cannot  be  ascertained  whether  a precipitation 
event  was  snow  or  rain.  Thus  the  model  in  Opus  is  only  approx- 
imate . 


Accumulation  of  Snow 

In  Opus,  snow  is  assumed  to  be  the  precipitation  form  when  the 
(typically)  regenerated  (or  measured)  sequence  of  daily  maximum 
and  minimum  temperatures  indicates  that  the  mean  daily  tempera- 
ture at  the  soil  surface  is  less  than  0 C.  Snowpack  is  accumu- 
lated from  any  precipitation  occurring  as  long  as  mean  daily 
temperature  remains  below  0 C. 
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Simulation  of  Snowmelt 


Snow  melt  is  calculated  on  days  when  a snowpack  exists  and  when 
daily  maximum  temperature  rises  above  0.  Figure  26  illustrates 
the  processes  that  Opus  considers  in  estimating  the  melt  and 
evaporation  of  snow.  As  with  rainfall,  a certain  amount  of 
snow  is  intercepted  on  any  plant  or  standing  dry  material. 

This  intercepted  snow  is  subject  to  early  evaporation  by  avail- 
able solar  and  thermal  energy.  Snowpack  evaporation  is  calcu- 
lated to  be  analogous  to  bare-wet-soil  evaporation,  with  PET 
obtained  using  the  estimated  albedo  for  snow-covered  surfaces. 

The  daily  effective  melt  temperature  is  apportioned  for  that 
part  of  the  day  that  is  estimated  to  have  temperature  above 
0 C,  assuming  sinusoidal  temperature  variation.  Melt  is 
assumed  to  come  from  two  potential  heat  sources:  air  heat 
convection  and  soil  diffusive  heat  flux.  Air  heat  melt  is 
estimated  using  a degree-day  method  modified  by  shade  and  snow 
depth  (Linsley  et  al . 1958).  Snow  depth  is  used  as  a simple 
surrogate  for  snow-heat  storage  and  ripening,  which  cannot  be 
properly  simulated  without  more  detailed  temperature  informa- 
tion and  heat- transfer  simulation  within  the  snow  cover.  The 
Opus  approximation  of  the  effect  of  snowpack' s latent  heat 
storage  is  illustrated  in  figure  27. 

Melt  due  to  soil-heat  flux  is  estimated  whenever  the  surface- 
soil  temperature  is  greater  than  the  assumed  snow  temperature 
(0  C) . Transfer  of  heat  from  soil  to  snow  is  assumed  to  be 
reduced  by  the  presence  of  surface  residue  and  mulch.  The 
combined  equation  for  estimating  daily  snowmelt  potential  (x^) 
(mm)  is 

xM  = 4. 5(1-0. 4F  2 )T  At  c.  + 10 . 8T  2k  (1-F  2)At  [106] 

M p m m f s*  s m 

in  which 

F is  relative  cover  shading  in  snow, 

T^  is  effective  mean  melting  temperature  (>0  C) , 

A?  is  effective  time  of  melting  (days) , 

cfm  is  exp( - .005Ws)2 , 

W is  snowpack  water  equivalence  (mm) , 

T 2 i-s  soil  temperature  at  10  mm  depth, 
k is  surface -soil  thermal  conductivity 
(mcal/cm/sec/  C) , 

F is  relative  mulch  cover  on  soil,  and 
A?  is  time  interval  (days) . 

The  first  term  on  the  right  of  equation  [106]  is  the  degree-day 
estimate  of  melt  from  atmospheric  heat.  The  second  term  is  the 
estimate  of  melt  due  to  thermal  flux  from  the  soil.  The  coef- 
ficient 4.5  is  the  degree-day  melt  coefficient  (Linsley  et  al . 
1958),  and  the  10.8  is  a unit  conversion  factor  for  the  latent 
heat  melt  for  heat  from  the  soil.  The  term  c^  (illustrated  in 
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fig.  27)  is  the  surrogate  effect  of  the  greater  coldness  of  the 
deeper  snowpacks.  ji  ^ 

Estimated  snowmelt  water  is  added  as  surface  input  during  cal- 
culation of  the  daily  soil  water  flow,  with  sinusoidal  distri- 
bution through  the  day  for  days  of  melt.  The  amount  that  is 
estimated  to  be  excess,  when  melt  rates  are  high,  is  treated  as 
runoff  in  the  surface -hydrology  section  of  the  model.  When 
rain  falls  on  a day  when  snow  is  melting,  the  two  inputs  are  t 

treated  in  sum  in  the  hydrologic  simulation.  z 

: 

* 


Figure  18. 

Two  examples  of  curve  number  relation  compared  with 
actual  rainfall-runoff  data.  Measured  daily  rainfall  and  runoff, 
when  plotted  in  this  way,  always  show  considerable  scatter,' 
which  the  more  simple  or  lumped  runoff  relations  cannot 
simulate  accurately. 
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PRECIPITATION,  mm 

Figure  19. 

The  SCS  method  for  treating  initial  soil-water  condition 
results  in  only  a very  small  relative  variation  in  predicted 
runoff.  The  two  curves  shown  are  the  CN  relations  for 
antecedent  conditions  I and  ll.The  case  simulated  is 
Watkinsville,  GA. 


The  Opus  breadpoint  hydrology  relates  infiltration  capacity 
(f)  to  infiltrated  depth  (I).  At  short  times  this  relation  is  an 
inverse  one  as  shown.  Actual  infiltration  rate  is  the  smaller  of 
f or  rainfall  rate  (r).  At  larger  times  or  at  larger  values  of  I,  the 
relation  becomes  asymptotic  to  the  constant  K . 


C 


Figure  20. 

The  probabilistic  method  in  Opus  assumes  that  the  Curve  Number  method  estimates  an  expected 
value  of  runoff.  At  any  level  of  rain  (P)  in  left  figure  (a,  b,  or  c)  there  is  actually  a random  distribu- 
tion of  runoff  (Q).  A two-parameter  beta  distribution  can  represent  a reasonable  distribution  for 
relative  runoff  amount,  as  shown  in  the  right  figure. 
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SOIL  MATRIC  POTENTIAL,  ^ 


The  effect  of  a surface  soil  crust  on  potential  gradient  in 
infiltrating  soil.  This  effect  is  important  in  simulating  the 


infiltration  relation  for  the  crusted  soil  case,  which  requires 
knowing  the  ultimate  value  of  f. 
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Figure  24. 

Two  curves  representing  hypotheses  used  by  Opus  to 
estimate  the  relation  between  soil  clay  content,  accumulated 
rainfall  energy,  and  formation  of  restricted  crust  layers  at  the 
surface  of  cultivated  soil.  1C  is  value  of  Kg  in  surface  layer  at 
any  time. 
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INFILTRATED  DEPTH,  I,  cm 

Figure  23. 

An  analytic  weighting  method  used  in  Opus  for  estimating 
infiltration  in  crusted  soils.  It  is  compared  here  with  the 
numerical  solution  of  Richards'  equation  for  the  same  case. 


NEW  SECTION  FOR 
DEPOSITION  CONDITION 


EVENT  HIGH  WATER 
AT  SECTION 


NEW  SECTION  FOR 
EROSION  CONDITION 


CALCULATED  STABLE  SECTION 
WIDTH  FOR  EVENT  HIGH  FLOW 


Figure  25. 

Section  changes  estimated  by  Opus  based  on  an  estimate  of 
altered  section  bottom  width  for  distribution  of  bottom  shear. 
Erosion  and  sedimentation  occur  in  different  amounts  at 
different  locations  along  a channel  or  furrow,  and  cause 
changes  in  the  furrow  geometry. 
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5.0 


Figure  26. 

Various  processes  that  affect  snowmelt  and  snow  evapora- 
tion. Snowmelt  is  assumed  to  include  loss  of  intercepted 
snow  from  vegetation.  Surface  accumulations  are  melted 
from  convective,  radiative,  and  soil  heat  energy  sources. 


Figure  27. 

Opus  approximation  of  effects  of  latent  heat  and  shading  on 
snowmelt. 
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6.  MICROBIOLOGICAL  AND  CHEMICAL  PROCESSES  IN  SOIL 


Transformations 

Nutrients 


of  Because  nutrients  are  a major  contributor  to  agricultural  non- 
point-source pollution,  it  is  important  to  accurately  simulate 
cycles  of  nutrients  in  the  agricultural  system.  These  cycles 
describe  the  exchange  of  nutrients  (including  those  in  applied 
chemicals  and  manures)  between  plants,  soil,  and  microbes. 

Along  with  the  transformation  cycles  of  phosphorus,  carbon,  and 
nitrogen,  the  movement  of  all  nutrients  in  tillage,  mixing,  and 
transport  in  soil-water  movement  should  be  considered  in  the 
nutrient  model.  All  of  these  processes  together  determine  the 
amount  of  movable  material  present  at  the  surface  when  runoff 
occurs,  and  the  amount  transportable  with  moving  soil  water. 

Organic  matter  is  central  to  the  cycling  of  plant  nutrients, 
and  the  simulation  of  nutrient  transformations  in  Opus  is  done 
by  a version  of  the  "Century"  organic -matter  model  of  Parton  et 
al . (1988b).  The  model  simulates  both  labile  and  stabilized 
organic -matter  fractions,  and  thus  simulates  the  nutrient- 
supplying  capacity  of  the  soil  organic  matter  (SOM)  as  well  as 
the  local  amount  of  mobile  nutrients.  A more  complete  discus- 
sion of  the  modeling  approach  and  the  background  for  the  model 
can  be  obtained  in  Parton  et  al.  (1988b). 

The  nutrient  model  contains  three  SOM  fractions  (fig.  28) , as 
follows : 

(a)  An  active  SOM  fraction  of  soil  carbon  (C) , nitrogen  (N) , 
and  phosphorus  (P) . This  consists  of  live  microbial  matter  and 
microbial  products,  along  with  partially  humified  soil  organic 
matter  with  a short  turnover  time  (1  to  5 y) . 

(b)  A slow  SOM  fraction  of  C,  N,  and  P that  is  physically 
protected  or  more  biologically  resistant  to  decomposition. 

This  fraction  has  an  intermediate  turnover  time  (20  to  40  y) . 

(c)  A passive  SOM  fraction  that  is  chemically  resistant  and  may 
also  be  physically  protected.  This  fraction  has  the  longest 
turnover  time  (200  to  1500  y) . 

Plant  residue  (above  and  below  ground)  is  divided  into  two 
fractions:  Structural  carbon  and  nutrient  pools  that  have  1- 

to  5-y  turnover  times,  and  metabolic  carbon  and  nutrient  pools 
that  have  0.1-  to  1-y  turnover  times  before  transfer  into  SOM 
pools.  Decomposition  products  from  these  pools  are  added  to 
the  SOM  pools  described  above. 

Decomposition  of  each  of  the  organic -matter  fractions  is  calcu- 
lated using  the  following  equation: 


d(m  ) . 
rc  l 

dt 


w 


<mrC>i 


[107] 
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where 


carbon  In  the  respective  pools  (gm/ra2), 
pool  index:  1 for  structural  and  2 for 
metabolic  soil  surface  litter;  3 for  structural 
and  4 for  metabolic  soil  litter;  and  5,  6,  and  7 
for  active,  slow,  and  passive  SOM,  respectively; 
time  (days) ; 

parameter  of  maximum  daily  decomposition  rate 
for  the  ith  state  variable; 
effect  of  relative  soil-water  content  on 
decomposition  (see  fig.  29  top);  and 
relative  effect  of  daily  mean  soil  tempera- 
ture on  decomposition  (see  fig.  29  bottom) . 


The  factors  of  temperature  and  water  content  were  developed 
based  on  laboratory  and  field  observations  (Parton  et  al . 
1987) . 


The  decomposition-rate  parameter  (k.)  is  constant  except  for  kx 
and  k3 , which  are  functions  of  the  lignin  content  of  the 
structural  material,  and  K6  (the  active  SOM  decay  rate),  which 
is  a function  of  the  soil  texture.  Decay  rates  of  surface  lit- 
ter are  slower  than  those  of  soil  litter,  based  on  the  assump- 
tion that  soil-moisture  conditions  are  consistently  less 
advantageous  for  decomposition  at  the  surface. 


The  model  assumes  that  all  pool  transfers  from  carbon  decompo- 
sition are  a result  of  microbial  activity  and  that  microbial 
respiration  is  associated  with  each  of  these  transfers.  The 
relative  amounts  of  loss  of  m ^ from  respiration  are  shown  in 
figure  28.  The  split  of  plant  residue  into  structural  and 
metabolic  fractions  as  a function  of  the  ratio  of  lignin  (m  ^) 
to  nitrogen  (m^)  was  determined  experimentally  by  incubation 
studies  (Pinck  et  al . 1950).  This  split  is  determined  using 
the  following  equation: 


f = 0.83  - 0.018 (m  T /m  M) 
rm  rL  rN 


[108] 


where  f is  the  fraction  of  residue  that  is  metabolic.  The 
fraction  that  is  structural  (f  ) equals  1 - f . The  split 
occurs  when  plant  residue  is  transferred  into  litter  material 
or  plowed  into  soil  residue.  The  fraction  of  structural  mate- 
rial in  the  litter  that  is  lignin  (f  ^)  affects  the  decomposi- 
tion rate  coefficients  kx  and  k3  in  the  following  manner: 


where  k, 
lc 


klc  exP(-3frL> 
k3c  exP(-3frL> 

and  k^c  are  appropriate  constants. 


[109] 

[HO] 
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These  relations  cause  the  decay  rates  of  structural  material  to 
decrease  as  the  relative  lignin  content  increases.  The  rela- 
tions are  based  on  the  assumption  that  the  microbes  can  more 
easily  decompose  certain  substrates  when  lignin  contents  are 
lower  (Melillo  et  al . 1984).  Values  of  parameters  in  equations 
[109]  and  [110]  were  obtained  from  laboratory  incubation  data. 

Laboratory  incubation  data  were  also  used  to  determine  func- 
tional forms  for  the  effect  of  soil  texture  (Sorenson  1981). 
Texture,  in  terms  of  sand  fraction  (fg)»  is  assumed  to  affect 
respiration  efficiency  (e  ) of  the  active  SOM  pool  and  its 
decay  rate  (k5),  as  follows: 

kc  = kc  (0.25  + 0.75f  ) [111] 

5 5c  s 

e = 0.17  + 0 . 68f  [112] 

r S 

The  relations  shown  above  have  been  verified  by  several  studies 
(Parton  et  al . 1987).  Determination  of  a value  for  k7  could 
not  be  done  by  laboratory  exercises  but  was  found  by  observa- 
tions on  stable  SOM  levels  at  several  sites  and  model- tuning 
procedures  (Parton  et  al . 1987). 

The  model  operates  on  a daily  time  step,  which  is  a modifica- 
tion over  other  versions  of  the  same  model.  Accumulations  of 
added  residue  material  and  losses  from  crop  harvesting  and 
grazing  are  accounted  for  each  day,  if  generated  by  the  manage- 
ment component.  The  location  of  the  SOM  pools  is  the  top  200 
mm  of  soil,  assumed  to  be  the  active  decay  region  of  soil. 
Usually  little  residue  decay  takes  place  below  this  level.  The 
mineral  N pool  is  subject  to  transport  through  this  zone  as 
water  moves  through  it.  This  is  accomplished  by  the  subsur- 
face-water-transport component  of  the  model  (ch.  4).  The  model 
is  linked  to  the  plant-growth  component  for  uptake  of  mineral  N 
and  P and  for  input  to  litter  of  residue  after  senescence,  or 
input  of  soil  residue  from  plowing. 


Submodel  for  Nitrogen 

The  nitrogen  model  has  the  same  flow  or  state  - transfer  struc- 
ture as  the  carbon  model  (see  fig.  30) . It  is  assumed  that 
most  N is  bonded  to  carbon.  Ratios  of  C/N  for  structural, 
active,  slow,  and  passive  soil  fractions  are  assumed  to  remain 
fixed  for  a given  site.  Ratios  of  C/N  for  structural  and 
active  SOM  are  set  at  150  and  8,  respectively,  and  the  C/N 
ratios  for  slow  and  passive  SOM  are  11.  These  ratios  are  based 
on  soil  pedon  data  analysis  and  other  data  (Parton  et  al . 

1987) . The  N content  of  the  metabolic  pool  is  allowed  to  vary; 
any  incoming  plant  N material  not  needed  for  the  structural  C/N 
ratio  goes  to  metabolic  N.  N transfers  are  stoichiometrically 
related  to  carbon  transfers,  with  N transfer  being  equal  to 


carbon  transfer  divided  by  the  C/N  ratio.  Transfers  between 
pools  can  result  in  mineralization  or  immobilization,  depending; 
on  the  C/N  ratios  of  the  donating  and  receiving  pools  and  also 
the  fraction  of  carbon  lost  as  C02  respiration  (30%  to  80%). 

Nitrogen  is  added  to  the  system  as  fertilizer  mineral  N,  as 
plant  residue  N (discussed  above) , or  by  biological  fixation. 
Fertilizer  N comes  in  bulk  as  an  addition  to  the  mineral  N 
content  (N03  or  NH4 ) of  a given  soil  layer,  usually  the  top  of 
the  soil. 

Other  N- transforming  processes  include  fixation  by  air  and 
plant  processes  and  also  nitrification  and  denitrification 
(Parton  et  al . 1988a).  Daily  fixation  from  the  atmosphere 
( itl,  ) (g/m2)  is  estimated  by 


"W  - CA  + VP)  l113' 

in  which  c^  and  b^  are  fixed  coefficients  and  P is  daily  pre- 
cipitation. Fixation  by  nitrogen- fixing  plants  is  assumed  to 
occur  to  the  extent  necessary  to  prevent  N stress  on  such 
plants . 


Nitrification,  which  changes  NH4  N to  N03  N,  is  assumed  to 
occur  continuously  in  response  to  soil  water  content,  tempera- 
ture, and  concentration  of  NH4 , as  follows: 

mNi  “ CNi  fW  fT  fNH  ^114^ 

in  which  the  f coefficients  are  factors  of  water,  temperature, 
and  NH4  concentrations,  respectively.  The  quantity  c^  is  a 
coefficient  with  a value  of  0.1  for  estimating  m^  in 
gm/m2/day.  Figure  31  illustrates  the  f functions  used  in 
equation  [114].  Denitrification  is  modeled  similarly,  with  an 
additional  soil -texture  factor  based  on  content  of  sand. 
Denitrification  causes  mineral  N03  to  be  volatilized  in 
oxygen-poor  situations.  Oxygen  poverty  is  usually  associated 
with  high  water  contents,  and  denitrification  is  assumed  to  be 
approximated  as  follows  (Parton  et  al . 1988b): 


mDi  Cdn  fW  fT  fSN  N°3 


[115] 


with  factors  f similar  to  those  for  nitrification  but  based  on 
nitrate  content  (N03)  (in  parts  per  million)  rather  than  on  NH4 
concentration.  The  quantity  f M is  an  approximator  of  texture 
and  is  defined  as  4.66  - 4.8fg  (f<,  = sand  fraction).  The 
coefficient  c^n  is  0.0004.  Figure  32  shows  the  factor  f^  for 
this  process  as  a function  of  water  content.  The  factor  f^  is 
the  same  as  for  nitrification.  Units  of  nu.  are  again  in 
gm/m2/day,  with  the  computed  amount  being  taken  from  the  soil's 
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Transformations 

Pesticides 


total  N03  content  and  lost  as  gas.  Amount  of  N volatilized  is 
assumed  to  be  5 percent  of  the  total  N mineralization  flows 
that  accompany  decomposition  of  metabolic  residue  and  SOM 
fractions . 

Submodel  for  Phosphorus 

The  model  for  phosphorus  is  similar  to  that  of  nitrogen,  with 
some  small  differences  (fig.  33)  (Parton  et  al . 1988b).  Ratios 
of  C/P  are  used  to  trace  P transformations  that  accompany  C 
transformations  between  the  various  pools.  The  primary  miner- 
al source  of  P is  weathering  of  soil  apatite.  During  weather- 
ing, labile  P is  taken  up  by  organisms  at  the  same  time  that 
secondary  and  occluded  forms  of  P are  produced  (fig.  33).  The 
SOM  pools  are  assumed  to  have  characteristic  C/P  ratios,  and 
the  amount  of  P moving  between  pools  is  determined,  as  for  N 
flows,  by  use  of  C/P  ratios.  Also  as  for  N,  P involved  in 
respiration  is  assumed  to  be  mineralized,  and  decomposition  of 
structural  residue  (which  is  low  in  N and  P)  is  assumed  to 
involve  immobilization. 


of  As  outlined  in  chapter  8,  pesticides  are  applied  in  several 
ways,  and  they  enter  the  simulation  by  deposition  on  plants, 
deposition  on  surface  soil  and  residue,  or  incorporation  into 
soil  at  some  depth.  The  portion  of  a pesticide  on  a plant  sur- 
face is  subject  to  both  environmental  decay  and  partial  washoff 
during  subsequent  rain.  Description  of  a pesticide  in  input 
data  includes  the  fraction  of  plant-adsorbed  pesticide  subject 
to  washoff  and  also  the  first-order  decay-rate  coefficient  for 
this  plant  fraction. 


Pesticides  within  the  soil  are  assumed  to  be  subject  to  envi- 
ronmental decay.  The  decay  rate  is  modified  by  temperature  and 
water  content,  but  is  assumed  to  be  described  by  a first-order 
rate  equation  (exponential  decay)  with  coefficient  kr  , as 
follows : 


ps 


dm  /dt  = k m [116] 
P'  ps  P 1 J 

in  which  nip  is  local  pesticide  mass.  The  effects  of  water  and 
temperature  are  assumed  to  be  described  by  the  Arrhenius  equa- 
tion (Walker  1974),  as  follows: 


kps  - kpsb  f„<*>  exp(bar/1.99(l/TKb-l/TKl)  [117] 


in  which 
T 

tK 

fKb 

w 


is  current  Kelvin  temperature, 
is  a reference  Kelvin  temperature, 
is  a water-activity  function  used  for  microbial 
degradation  (see  fig.  29), 
is  a constant,  and 

is  decay  coefficient  at  temperature  . 
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The  coefficient  b is  calibrated  if  data  are  given  for  a mea- 
sured K k at  a given  temperature  or  water  content,  or  else  is 
assumed^lo  have  a default  value  taken  from  published  values 
(Walker  1974) . Without  better  and  more  widely  available 
knowledge  of  the  individual  processes,  Opus  does  not  attempt  to 
separate  the  various  processes  of  environmental  degradation  and 
transformation  of  pesticides  (such  as  hydrolysis  and  photochem- 
ical action) . Decay  is  assumed  to  be  highly  related  to  micro- 
biological activity,  so  it  follows  that  a pesticide  leached 
beyond  the  depth  of  active  decay  (approximately  200  mm)  is 
decayed  at  a fraction  of  the  decay  rate  of  pesticides  in  the 
active  region. 


Pickup  of  Chemicals 
by  Rainfall  and 
Runoff 


Several  processes  are  simulated  in  the  Opus  model  in  connection 
with  the  production  of  chemicals  in  runoff  water.  The  simula- 
tion methods  are  designed  to  operate  equally  well  with  both  of 
the  sediment  transport/hydrology  options  available.  This  dic- 
tates that  the  methods  are  lumped.  It  is  doubtful  that 
detailed,  spatially  distributed  calculations  have  general  value 
because  of  many  nonuniformities  and  our  imprecise  knowledge  of 
the  exact  mechanisms  in  the  various  processes. 

The  following  paragraphs  outline  the  methods  used  in  Opus  to 
estimate  the  load  of  nitrate,  phosphorus,  and  ammonia  in  runoff 
water,  and  the  amount  of  each  that  washes  into  the  soil  from 
the  surface.  These  constituents  are  treated  identically,  with 
few  noted  exceptions.  The  following  discussion,  for  simplici- 
ty, uses  labile  phosphorus  (P)  as  an  example. 


The  movement  of  both  soil  and  surface  water  is  involved  in  the 
transport  of  chemicals  in  the  agricultural  environment.  How- 
ever, the  occurrence  of  a runoff -producing  rain  has  the  primary 
role  in  actually  moving  chemicals  off  the  field.  Nonrunoff 
rainfall  events  can  also  move  chemicals  from  plant  and  residue 
surfaces  into  the  soil  surface. 


Leaching  of  Canopy  Chemicals 

For  all  rainfalls  that  exceed  the  interception  storage  of  the 
current  canopy,  a certain  small  fraction  of  the  nitrate  and 
phosphorus  contained  in  the  plant  is  assumed  to  wash  off  the 
plant  with  the  rainwater.  For  living  plants,  this  fraction  is 
assumed  to  be  0.02;  for  standing  residue,  or  plants  undergoing 
senescence,  it  is  assumed  to  be  0.03  (Schrieber  1990).  If  the 
period  since  the  last  rainfall  or  washing  event  [At^  (days)]  is 
greater  than  1.0,  the  fraction  is  reduced  by  a factor  (f  ) as 
follows: 

f =1  -exp(6At.)  [118] 

pw  a 
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In  addition,  the  amount  of  plant  leachout  is  proportional  to 
the  percent  cover  that  the  crop  represents  compared  to  the 
field  surface  (F  ).  Thus,  the  estimated  field  value  of  labile 
phosphorus  leachld  from  the  canopy  [mp(cl) , in  kg/ha]  is 

mp(cl)  = 0.02mp(c)fpwFp[LFL/FLT] 

in  which 

mp(c)  is  P content  in  crop  (c)  (kg/ha), 

is  leaf  area  index  for  that  crop,  and 

FT  is  total  field  Leaf  Area  Index  (F  ) for  all  plants. 

LT  -L* 

Leaching  of  Residue 


[119] 


Rainwater  also  leaches  decomposition-product  nutrients  from 
residue  on  the  surface  of  the  soil.  This  process  is  assumed  to 
be  similar  to  the  clarification  of  a simple  reservoir.  Adding 
water  to  a solution  in  an  overflowing  reservoir  results  in 
exponential  decay  of  concentration,  assuming  a fixed  mixing 
efficiency.  The  leaching  constant  for  residue  is  assumed  to  be 
.05  per  nun  of  rainfall.  The  amount  of  P in  the  residue  that  is 
leached  in  a given  storm  [m  (rsl) ] is  assumed  ( ye  Schrieber 
1990)  to  be 

mp(rsl)  - mp(rs)(l-exp[- .05V  ] ) [120] 

where 


mp(rs) 

V 

P 


is  current  content  of  P in  surface  residue 
(kg/ha) , and 

is  rainfall  passing  through  residue  (mm) . 


i 


Routing  Chemicals  by  Infiltration  and  Runoff 

Phosphorus  or  other  nutrients  from  the  above  sources,  plus 
natural  nitrate  in  rainwater,  will  either  enter  the  soil  or 
leave  the  catchment  in  runoff.  The  path  of  a solute  is  esti- 
mated according  to  the  division  of  rainfall  in  the  infiltration 
process.  The  first  part  of  the  rainfall  will  all  infiltrate 
into  the  soil.  This  portion  (I  ) is  then  distributed  downward 
from  the  soil  surface  (along  w?th  chemicals  contained  in  the 
water)  to  the  depth  reached  by  the  wetting  front  at  ponding 
(z  ) . A more  detailed,  spatially  distributed,  and  dynamic 
st&dy  has  been  made  by  Havis  (1986).  The  model  here  is 
necessarily  lumped. 
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Rainfall  and  washed  solutes  reaching  the  soil  surface  after 
ponding  will  be  divided  between  infiltration  and  runoff,  but 
all  incoming  rain  will  interact  with  the  flowing  surface  water. 


First,  the  concentration  of  chemicals  in  the  surface  active 
layer  is  calculated  at  the  time  when  runoff  begins,  based  on 
the  amount  I , which  has  washed  through  it.  Then  the  rainwater 
that  exceeds^I  interacts  with  this  surface  zone  as  a simple 
reservoir  duriRg  runoff.  For  chemicals  such  as  phosphorus  and 
ammonium,  this  interaction  assumes  the  behavior  of  an  equilib- 
rium adsorption  isotherm.  The  depth  of  the  interactive  soil  is 
related  to  soil  disturbance,  porosity,  and  infiltration  charac- 
teristics. The  depth  is  user  specified  as  a fraction  (f  ) of 
the  10-mm-deep  surface  layer. 


The  mixing  interaction  is  assumed  to  be  similar  to  that  in 
soil-water  transport.  The  estimated  concentration  of  runoff 
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122] 


in  which 


M = total  water  for  equivalent  dissolved  labile 
P at  equilibrium  in  active  layer  (mm) , 

= 105(1)  + 1000  Kd  M(l) , 

K^p  = equilibrium  adsorption  isotherm  for  P (liter/kg), 

M(l)  = mass  of  soil  in  top  10  mm  (kg/ha), 

5(1)  = water  content  in  first  (10-mm)  layer  after 

infiltration, 

r = rainfall  reaching  soil  after  surface  ponding 

(mm)  . 


Finally,  the  amount  of  P in  this  example  is  assumed  to  equili- 
brate with  the  transported  sediment  according  to  the  equilibri- 
um adsorptivity  described  by  K^p,  which  partitions  the  P in 
runoff  between  that  adsorbed  on  sediment  particles  and  that  in 
solution . 


Pesticides  enter  runoff  water  by  very  similar  assumed  mecha- 
nisms. Some  pesticides  may  contribute  part  of  their  concentra- 
tion by  washoff  from  plants.  In  Opus,  this  portion  is  kept 
separate  from  interactions  between  runoff  water  and  soil 
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pesticides,  for  simplicity.  An  additional  complexity  for 
pesticides  enters  when  kinetic  adsorption  dynamics  is  assumed. 
For  that  reason,  the  equilibrium  distribution  used  in  equations 
[121]  and  [122]  cannot  be  treated.  An  estimate  is  first  made 
of  the  average  concentration  of  the  dissolved  and  adsorbed 
amounts  of  each  pesticide  in  the  upper  layer  during  the  runoff 
time.  A weighted  average  is  used.  The  weighting  is  propor- 
tional to  the  ratio  of  infiltrated  water  to  the  water  content 
of  the  first  layer.  More  infiltrated  water  will  more  heavily 
weight  the  concentrations  at  the  end  of  the  runoff  event.  Then 
the  pesticide  is  assumed  to  interact  with  the  runoff  water, 
using  a surface  active  fraction  as  above,  and  assuming  a 
reduced  equilibration  mass  fraction  based  on  the  kinetic 
adsorption  rate  factor  (u)  and  length  of  rainfall  (t) , as 
follows : 


f^  = exp  [ t(  . 01-j/)  ] 


[123] 


The  factor  f^  is  multiplied  by  the  surface  interactive  fraction 

(f  ) to  obtain  the  equivalent  amount  of  soil  from  which  the 

adsorbed  pesticide  will  equilibrate  with  the  runoff  water.  For 

the  equilibrium  option,  f^  is  1.  The  total  pesticide  from  the 

soil  available  for  runoff  interaction  (M  ) is  then 

pt 


M 


pt 


- W + V,”*1) 


[124] 


where 


M 

ps 


is  the  mass  of  pesticide  in  soil  solution  in  the 
upper  layer  (1) , and 


is  adsorbed  concentration  of  pesticide. 

Writing  an  equation  for  this  amount  equilibrated  with  the  total 

runoff  water  depth  Q (mm) , and  solving  for  the  concentration  of 

runoff  dissolved  pesticide  (C  ) , we  obtain 
r ro 

C = M /[l/cu  + f fuM(l)KJ  [125] 

ro  pt  L ' h z b dJ  L J 

in  which  c^  is  the  conversion  factor  defined  below  equation 
[37].  The  amount  of  pesticide  dissolved  in  runoff  in  kg/ha  is 
this  concentration  times  Q divided  by  the  conversion  factor  c^. 
The  concentration  on  sediment  (C  ) results  from  surface  soil 
that  has  released  pesticide  into^runoff,  according  to  fraction 
f^,  plus  transported  particles  with  adsorbed  material,  as 
follows : 


C 

pm 


= K ,C  f, 
d ro  b 


+ (l  - fb)c 


a 


[126] 
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where  C refers  to  effective  surface-soil-adsorbed  concentra- 
tion  during  the  runoff  event,  as  discussed  above.  This  concen- 
tration is  then  multiplied  by  the  enrichment  ratio  and  the 
sediment  mass  to  obtain  the  estimated  amount  of  pesticide  of 
each  species  that  leaves  the  catchment  on  sediment.  The 
amounts  adsorbed  and  dissolved  in  surface  soil  are  then  adjust- 
ed for  this  loss  by  mass  balance. 


Figure  28. 

Diagram  for  the  century  model  of  the  organic  carbon  cycle 
(used  in  Opus)  and  the  flow  of  carbon  between  the  various 
pools.  (Reproduced  with  permission  from  Parton  et  al. 
1988b.) 
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Figure  29. 

Relative  effect  of  soil  moisture  (top)  and  soil  temperature 
(bottom)  on  microbial  decay  process  for  the  Century 
model  for  soil  carbon  and  nutrient  cycles. 
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Figure  30. 

Diagram  for  the  Century  model  for  organic  nitrogen  cycle 
and  the  flow  paths  for  exchange  of  N between  pools. 
(Modified  from  Parton  et  al.  1988b.) 
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Figure  31. 

Functions  for  the  effects  of  the  soil  water  content  (upper), 
soil  temperature  (middle),  and  soil  ammonia  (lower)  on 
relative  rates  of  nitrification. 
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Figure  32. 

Relative  effect  of  water  content  on  rate  of  denitrification  in 
the  Century  model  for  nitrogen  cycling. 


Temperature  Factor 


Figure  33. 

Diagram  for  flows  between  various  pools  in  the  Century 
model  for  phosphorus  in  Opus.  (Reproduced  with  permission 
from  Parton  et  al.  1988b.) 
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7.  SIMULATION  OF  PLANT  GROWTH 


Opus  uses  a mechanistic  plant  model  to  simulate  the  growth  of 
plants  or  crops  in  response  to  four  major  factors  related  to 
plant  development:  radiation,  nutrients,  temperature,  and 
water  availability.  The  model  may  be  used  for  either  annual  or 
perennial  crops,  and  it  responds  to  grazing,  harvest,  and 
stresses  in  a relatively  realistic  manner.  It  is  distantly 
related  to  the  plant-growth  model  in  the  EPIC  model  (Williams 
et  al.  1984),  and  many  of  the  Opus  concepts  were  inspired  by 
early  versions  of  the  EPIC  plant  model  (Williams,  personal 
communication) . 


Growth  Model 


In  general,  the  plant  model  can  be  expressed  as  a relation 
between  various  factors  and  the  incremental  production  of  plant 
material  (Apm)  in  a time  period  (At),  as  follows: 


Apm  =c  f f f r R.  At 
e e a m i 

in  which 


[127] 


c = photosynthetic  conversion  coefficient  for 
production  of  plant  material  from  radiant 
energy  (kg/ha/langley ) , 

fg  = factor  (0  to  1)  for  energy-conversion  efficiency, 
f - factor  (0  to  1)  for  plant  age, 

3. 

f = factor  (0  to  1)  for  plant  size, 

r = stress  factor  (0  to  1),  considering  water,  nutrients, 
and  temperature  stresses, 

= daily  radiation  (langleys) , 

At  = time  increment  (days) . 


The  value  c is  a time  - invariant , crop-related  parameter,  with 
the  conversion  efficiency  varying  over  a limited  range  from 
plant  to  plant.  The  values  f and  f are  functions  that  auto- 
regulate  the  plant  growth;  they  are  illustrated  in  figure  34. 
The  value  f slows  the  plant  growth  as  degree -day  age  approach- 
es maturity,  and  f slows  the  plant  growth  as  maximum  genetic 
size  is  approached.  These  functions  are  mechanistic  rather 
than  experimentally  determined,  but  they  produce  a plant 
response  that  imitates  observed  plant  behavior. 


Self-Dependent  Growth 

The  value  f reflects  change  over  time  in  the  plant's  ability 
to  intercept  radiation  and  photosynthesize  new  material  as  the 
leaf  area  increases  toward  that  at  leaf  maturity.  Growth 
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depends  on  the  amount  of  radiation  that  can  be  intercepted  by 
the  leaves  and  thus  on  the  amount  of  leaf  area.  The  relation 
of  leaf-area  index  (F^)  to  leaf  mass  is  discussed  below. 
Because  of  shading  and  the  variable  nature  of  maximum  (re- 
ferred to  here  as  FT.,),  the  relation  of  F to  FT  \FT . . is  direct 
at  small  values  of  F , but  is  asymptotic  to  the  upper  limit  of 
1.  The  equation,  illustrated  in  figure  35,  is 

fe  - 1 - exp(-3F h/Tm  - ve  - c.)  [128] 

In  perennial  plants  and  some  multiple -harvest  plants  (such  as 
alfalfa),  root  mass  can  be  large  relative  to  aboveground  mass, 
and  stored  energy  can  be  used  to  produce  leaves  and  plant 
material.  This  is  modeled  here  by  a factor  v (defined  below), 
which  boosts  growth  efficiency  above  that  limited  by  the  leaf 
area.  The  factor  c.  acts  similarly  to  represent  the  initial 
production  of  leaves  from  seeds.  This  factor  shrinks  to  zero 
as  the  plant  grows.  It  is  mathematically  necessary  to  initiate 
seedling  growth.  In  mathematical  terms,  equation  [127]  is  in 
simplest  form  a linear  differential  equation  at  small  values  of 
pm  (or  F^) , and  cannot  exhibit  growth  for  pm  = 0 , without  a 
value  for  c..  Equation  [128]  reduces  to  f = F + v + c . at 
small  values  of  F^.  6 

As  leaf-area  index  increases,  growth  rate  can  increase, 
expressing  a self-dependent  growth  rate.  The  parameter  on  the 
curves  in  figure  35  indicates  the  simulated  effect  of  excess 
root/leaf  ratio.  The  value  v is  the  ratio  of  mass  of  roots 
to  mass  of  leaves  in  excess  of  that  called  for  by  a normal 
root/leaf  ratio.  Positive  values  of  v can  occur  under  hay  or 
alfalfa  harvest  or  grazing.  This  ratio  also  expresses  an  early 
spring  start  for  perennials  with  an  established  root  system, 
compared  with  plants  starting  from  seed.  A value  of  vg  = 0 
indicates  normal  root/leaf  ratio  for  annuals. 


Growth-Limiting  Stresses 

The  stress  factor  (r)  is  the  minimum  of  individually  considered 
stress -coefficient  approximations  for  water,  temperature,  and 
nutrients.  It  is  1.0  for  no  stress  and  0 for  complete  cessa- 
tion of  plant  growth.  The  lower  the  value  of  this  factor,  the 
greater  the  stress  the  plant  is  experiencing  in  that  time 
period.  The  time  period  for  growth  computation  is  typically  1 
day . 
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Water  Stress 


The  water  stress  factor  (r  ) is  1 for  net  contents  of  root-zone 
water  down  to  20  percent  greater  than  15-bar  water  content 
(015),  and  declines  linearly  between  that  point  and  01S,  as 
follows : 

r =1.0 

r”  - (9  - 015)/( 0.281S) 

T = 0. 

w 

Nutrient  Uptake  and  Stress 

The  nutrient  stress  factor  (r^)  is  a function  of  the  ratio  of 
actual  to  potential  use  of  plant  nitrogen.  Phosphorous  stress 
is  not  simulated.  Plants  maintain  a given  N/P  ratio,  which  is 
used  in  the  residue-decay  model  (ch.  6).  Potential  use  of 
plant  nitrogen  is  obtained  from  information  on  the  N content  of 
the  plant  at  various  growth  sizes,  as  illustrated  in  figure  36. 
Plants  typically  contain  more  N at  emergence  than  at  maturity. 
The  relation  of  fractional  N content  to  relative  size  [n(d) ] 
can  be  expressed  as  a three -parameter  function,  based  on  N 
contents  at  emergence  and  maturity,  plus  a shape  factor,  as 
follows : 

n(d)  = n + (n  - n )exp( -c  ,d)  [130] 

m om  r d 

in  which 

d is  dimensionless  plant  size  expressed  as  ratio  of 

dry  matter  or  plant  material  (pm)  to  potential  maximum 
plant  material  (ppm) , = pm/ppm 

nQ  is  content  of  plant  nitrogen  at  emergence, 

n is  a parameter  set  so  that  n(d)  equals  mature 
N content  for  d = 1 , and 

c^  is  a coefficient  controlling  decay  rate  of  the 
function. 

This  allows  estimation  of  the  potential  content  of  plant  N at 
any  point.  If  root-zone  soil  N is  insufficient  to  meet  the 
daily  N demand  thus  estimated,  the  N stress  factor  (r^)  is 
calculated  as 


8 > 1.2 eXb 

8lb  < 8 < 1.2 8lb  [129] 

^ $ 1 6 


0 actual  N use 

j — 
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1 


[131] 


Allocation 

Material 


It  is  assumed  that  plant  roots  are  able  to  selectively  garner 
nitrogen  from  soil  water  if  the  N demand  by  the  plant  exceeds 
that  amount  present  in  the  daily  water  uptake.  Plant  use 
cannot  exceed  the  nitrate  and  ammonia  in  the  soil  water  ambient 
to  the  roots. 


Temperature  Stress 

The  plant- temperature  stress  factor  (r^)  is  a dimensionless 
function  of  the  daily  mean  air  temperature  (T) , the  minimum 
growth  temperature  (T^) , and  the  optimum  growth  temperature 
(T  ) (illustrated  in  fig.  37).  Several  types  of  curve  func- 
tioRs  are  candidates  to  represent  observed  plant- temperature 
responses,  and  plants  naturally  differ  somewhat.  The  func- 
tion used  in  Opus  may  not  exactly  represent  every  type  of 
plant.  The  normalized  temperature  (T^)  used  in  figure  37  is 
defined  on  the  abscissa. 


of  Plant  Plant  material  produced  by  the  model  thus  described  is  divided 
among  leaf/stalk,  root,  and  fruit/seed  material  according  to 
the  plant's  relative  age  and  size.  A general  division  function 
is  presently  used  for  all  plants,  and  is  illustrated  in  figure 
38.  Relative  root/leaf  ratio  is  a function  of  relative  mass 
(lower  scale) , and  relative  seed  or  fruit  production  is  a func- 
tion of  relative  age  in  degree-days  (upper  scale). 

The  plant  material  in  leaves  and  roots  must  be  interpreted  to 
estimate  leaf  area  (or  F^)  and  rooting  depth  to  complete  the 
plant  description.  The  relation  between  the  mass  of  plant 
material  in  leaves/stems  (pm^)  and  the  leaf  area  index  (F.^)  is 
assumed  to  follow  a nonlinear  function  as  follows: 

FL  ” FLM  sin[*/2  Pmiv/ppnilv^  t132! 

in  which 

is  maximum  potential  leaf-area  index  for  the  plant 
type , and 

ppm^v  is  maximum  potential  leaf/stem  dry  matter  for  the 
plant  type. 

The  plant  material  allocated  to  root  mass  (pm^)  results  in  a 
gradual  increase  in  depth  of  root  penetration  into  the  soil 
(Zr)  up  to  a maximum  potential  root  depth  (Z  ) , as  follows: 

Zr  = 10  + (Z  -10) |l-exp[ -2(exp[l . Spm^/ppm^] -1) ]|  [133] 
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Interaction  of 
Plant  Growth 
Factors 


in  which  ppm  is  potential  maximum  root  mass  for  the  plant. 
This  function  includes  the  allowance  for  planting  of  the  seed 
for  annual  crops  at  least  10  mm  deep. 


If  the  leaf - intercept  area  were  directly  proportional  to  plant 
dry  matter  and  if  the  age  and  senescence  factor  are  ignored, 
equation  [127]  would  be  a linear  differential  equation. 

Growth  could  not  initiate  without  the  presence  of  c.,  which  may 
be  thought  of  as  representing  the  ability  of  the  plant  to 
produce  leaves  from  seed  material.  As  described  above,  the 
plant  material  produced  is  divided  between  root  and  nonroot 
(leaf+stem,  and  fruit)  material.  If  we  assume  for  the  moment 
that  factors  f , f , and  R.  are  constants,  and  that  f is 
proportional  to  plus  the  value  c^  (as  is  true  at  small  F^) , 
the  equations  [13Z]  and  [133]  can  be  combined  and  simplified  to 
illustrate  the  underlying  basic  growth-rate  function  as 


dM 

dt 


B [ sin(M^)  + b] 


[134] 


in  which  (recalling  that  pm  is  plant  material  and  ppm  is 
potential  maximum  plant  material) 


is  n/2 (pm/ppm) , 

B is  n/2 (Y/ppm) , 

b is  c./Y,  and 

Y is  fLf  R.c  . 

a m l e 

The  solution  of  equation  [134]  has  the  form 


1 , b tan(m*/2)  -f  1-v  . 1+v 
B v n b tan(m^/2)  + 1+v  1-v 


in  which 


v = /1-b1 2 

Although  equation  [135]  includes  assumptions  neglecting  senes- 
cence factors,  it  can  help  obtain  an  index  of  the  optimum  time 
to  maturity  for  a set  of  plant  parameters,  by  solving  for  t 
when  reaches  its  peak  of  n/2.  Notice  that  although  the 
solution  is  scaled,  the  time  is  inversely  proportional  to  B and 
is  thus  proportional  to  ppm.  This  is  reasonable  because  a 
fixed  amount  of  radiant  energy  is  available,  which  causes  larg- 
er plants  to  reach  maturity  at  later  dates  than  small  ones,  all 
other  factors  being  equal. 
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Senescence 


After  maturity,  failing  a harvest  operation  that  removes  the 
plant,  active  leaves  will  turn  to  "standing  dry  matter"  at  a 
rate  specified  as  a parameter  by  the  user.  Roots  of  annual 
plants  are  converted  to  soil  residue  when  active  leaf  area 
falls  to  zero.  For  perennial  plants,  roots  are  assumed  to 
recycle  and  regenerate  10  percent  of  their  mass  over  winter. 

No  doubt  this  figure  actually  varies  widely  with  plant  type  and 
climate . 

In  summary,  the  major  parameters  needed  by  this  model  to 
describe  a plant  include  the  following: 

(a)  seasonal  lifetime  in  degree-days  to  onset  of  senescence, 

(b)  age  at  emergence  in  degree-days  (measured  from  planting  for 
annuals  or  from  January  1 for  perennials), 

(c)  potential  production  of  total  plant  material  for  unstressed 
condition  (r=l) , 

(d)  N content  with  growth  stage  (three  parameters), 

(e)  potential  weight  of  fruit/seeds  at  maturity, 

(f)  maximum  depth  of  rooting,  and 

(g)  conversion  efficiency  of  photosynthetic  energy. 


Discussion 

Like  some  other  parts  of  the  Opus  model,  the  plant  model  is 
well  behaved  (or  bounded  in  its  response  to  inputs)  but  also 
difficult  to  exhaustively  test  because  of  the  paucity  of  data 
on  the  growth  response  of  plants  to  known  amounts  of  stress. 

Few  data  exist  for  other  than  a few  crops  for  the  simultaneous 
values  of  plant  mass,  leaf  area,  water  content,  and  nitrogen 
available  in  the  root  zone.  Typically,  experimental  objectives 
are  limited  to  only  a few  of  these  factors.  The  model 
described  above  has  been  compared  with  various  types  of  plant 
data  from  several  climates  around  the  continental  United 
States,  and  parameters  are  given  in  the  User  Manual  (vol.  II). 
That  tabulation  should  guide  the  user  in  estimating  values  for 
related  crops  and  climate  zones. 


Example  of  Plant  Figure  39  illustrates  the  plant-model  simulation  of  a hypothet- 

Growth  Simulation  ical  corn  crop,  with  only  temperature  stress  imposed.  The 

combined  action  of  the  factors  f , f , and  f is  manifest  in 

e m a 
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the  shape  and  timing  of  the  growth  pattern  that  results.  The 
same  corn-parameter  set  produces  significantly  less  when  the 
stresses  of  limited  water  and  nutrients  are  included,  as  illus- 
trated in  figure  40.  This  particular  year,  1973,  was  a rela- 
tively dry  one.  Application  of  the  plant  model  to  a multiple- 
harvest-grass  field  crop  is  illustrated  in  figure  41. 


RELATIVE  AGE  (Deg -Days)  OR  RELATIVE  SIZE 


Figure  34. 

The  relative  growth  factors  f and  f , representing  relative 
aging  and  relative  size  used^n  the  mechanistic  plant  model. 


Figure  35. 

Illustration  of  equation  [128],  As  leaf  area  index  increases, 
the  relative  ability  of  the  plant  to  intercept  radiation  and 
produce  new  plant  material  increases  at  a decreasing  rate, 
with  an  upper  limit  of  1 . For  perennial  plants  or  cropped 
grasses,  energy  stored  in  roots  can  effectively  increase  this 
rate;  Vg  is  relative  proportion  of  roots  to  aboveground 
material  in  excess  of  normal  ratio. 


Figure  36. 

Concentration  of  plant  nitrogen  as  a function  of  relative  size 
or  growth  stage  for  several  common  crops.  Symbols  are 
from  published  data:  W = wheat  (Boatwright  and  Haas 
1961),  CN  = corn  (Hanaway  1962),  SB  = soybeans 
(Hanaway  and  Weber  1971),  SBF  = fertilized  soybeans 
(Hanaway  and  Weber  1971),  CT  = cotton  (Bassett  et 
al.  1 970). 


T*  = <T-Tb)/(Top-Tb> 

Figure  37. 

Relative  effect  of  temperature  on  growth  rate  of  plants, 
based  on  a scaled  or  relative  temperature.  Tb  is  lowest  and 
T is  optimum  temperature  for  plant  growth. 
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RELATIVE  PLANT  AGE 
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RELATIVE  PLANT  SIZE 

Figure  38. 

Division  of  plant  material  by  the  plant  growth  mode!  in  Opus. 
Plant  material  produced  by  the  plant  is  divided  between 
leaves,  roots,  and  fruit  according  to  the  relative  size  of  the 
plant  or  the  relative  age  of  the  plant  (upper  scale)  in  the  case 
of  fruit. 


Plant  model  simulation  of  ideal  growth  (no  stress)  of  a corn 
crop  in  Georgia. 


Figure  40. 

The  same  simulation  as  in  figure  39  except  that  water  stress 
was  included,  reducing  growth  and  yield. 


Figure  41. 

Example  Opus  plant  model  simulation  of  the  growth  of  a 
multiply  harvested  grass  crop.  This  figure  illustrates  the  rapid 
regrowth  feature  where  root  resources  are  available. 
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S.  SPECIFYING  MANAGEMENT  OPERATIONS 


Several  major  options  for  the  use  and  treatment  of  land  are 
specified  as  part  of  the  field  description,  discussed  in  chap- 
ter 2.  These  include  the  direction  of  tillage  and  the  use  of 
terracing,  an  impoundment,  and/or  various  grass  buffer  strips. 
Other  management  operations,  especially  those  with  specific 
dates  of  application,  are  discussed  below.  Management  is 
assumed  to  occur  on  a multiyear  rotation  basis. 

Management  operations  are  specified  by  a list  of  operations  and 
their  characteristics  plus  a schedule  of  operation  amounts  and 
target  dates.  The  list  describes  the  major  physical  changes 
that  each  operation  effects  in  the  condition  of  the  field  or 
crop.  The  schedule  specifies  the  date  on  which  that  operation 
should  occur  and  any  variable  amounts  involved.  Each  year  in 
the  rotation  cycle  has  a schedule,  and  the  rotation  period  can 
include  as  many  as  5 years. 


Lists  of  Choices  Opus  allows  users  to  choose  from  among  crops  to  grow,  tillage 

procedures  to  use,  and  pesticides  and  animal  wastes  to  apply. 
Certain  characteristics  of  members  of  these  four  lists  are 
specified  and  remain  fixed,  independent  of  when  the  farming 
operations  involving  the  choices  are  implemented. 


List  of  Crops 

The  set  of  parameters  representing  a crop  or  plant  (as  given  in 
ch.  7 in  connection  with  the  plant  model)  are  listed  for  each 
plant  species  that  is  grown  during  the  management  rotation 
period.  For  a perennial,  a mix  of  perennials,  or  an  unmanaged 
catchment,  the  rotation  period  is  1 year,  and  the  plant  charac- 
teristics are  listed  for  the  major  plants  or  plant  families 
found  on  the  catchment.  A maximum  of  four  plant  species  may  be 
treated  simultaneously.  When  an  annual  is  planted,  the  parame- 
ters in  this  list  are  used  by  the  plant  model  to  simulate  its 
growth.  The  crop  parameters  list  includes  specification  of  a 
type  for  each  plant.  The  plant  types  are  annual,  harvested 
perennial,  grazed  perennial,  and  annual  meadow  or  no -till. 

This  type  controls  how  the  plant  model  operates  and  how 
instructions  for  planting  and  harvesting  or  grazing  are  inter- 
preted. These  features  are  more  fully  explained  in  the  User's 
Manual  (vol.  II). 


List  of  Tillage  Operations 

Mechanical  operations  in  the  field  are  classed  into  four  types, 
numbered  as  follows : 
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(1)  Planting : A plant  seed  or  seedling  is  placed  in  the  ground 
at  a specified  depth,  during  which  some  mixing  of  the  soil  may 


occur.  A row  spacing  is  specified,  which  may  not  be  changed 
until  harvest.  Plant  degree-day  "timers"  (ch.  7)  are  started 
on  the  planting  day. 

(2)  Cultivation:  The  field  soil  is  disturbed  with  a specified 
mixing  efficiency  to  a specified  "plow"  depth.  Plant  material 
standing  at  the  time  is  assumed  to  be  unaffected.  Row  spacing 
is  unchanged,  but  furrow  depths  may  be  changed.  Any  shape 
changes  from  past  erosion  or  deposition  are  eliminated.  These 
characteristics  distinguish  cultivation  from  plowing  opera- 
tions. Mixing  includes  all  soil  constituents  down  to  the  given 
"plow"  depth. 

(3)  Harvesting : A specified  proportion  of  the  plant  is 
removed,  including  at  least  the  seed/fruit  portion  existing  at 
harvest  time.  The  soil  is  not  mixed,  and  row  spacing  and 
furrow  depth  are  not  affected.  A special  code  allows  harvest- 
ing of  the  plant  plus  roots  for  root  crops.  For  grazed  crops, 
this  operation  specifies  time  interval  (days)  and  rate  of 
removal  of  leaf  and  stem  material  by  animal  grazing. 

(4)  Plowing:  Characteristics  for  plowing  operations  are  simi- 
lar to  those  for  cultivation  except  that  live  or  dead  standing 
plant  material,  surface  litter,  and  plants  now  growing  are  also 
mixed  into  the  soil  to  a given  depth  with  a given  efficiency. 
Row  spacing  may  be  changed  arbitrarily. 

Data  for  operation  types  1 and  3 must  include  naming  a crop 
from  the  crop  list.  If  harvest-operation  type  is  specified  for 
a date  with  a crop  that  has  been  designated  a grazed  perennial, 
then  it  is  assumed  that  a period  of  grazing  is  started  on  that 
day  and  extends  a specified  number  of  days  with  a specified 
grazing  rate. 

Mixing  that  accompanies  management  types  1,  2,  and  4 causes 
major  changes  in  the  soil  material.  Soil  bulk  density  is 
returned  to  a minimum  value.  A mixing- efficiency  parameter 
determines  the  degree  of  mixing,  so  that  soil-water  content  and 
all  dissolved  and  adsorbed  materials  (nitrate,  labile  phos- 
phate, pesticides,  and  residue  material)  will  be  mixed  to  the 
specified  efficiency.  It  is  assumed  that  the  plowing  does  not 
extend  significantly  beyond  the  200-mm  depth  of  microbial 
activity  and  that  nutrient-model  pools  are  not  altered.  It  is 
also  assumed  that  plow  depth  does  not  exceed  the  depth  at  the 
upper  soil  horizon:  Opus  cannot  predict  the  results  of  mixing 
two  different  soil  types. 
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List  of  Pesticides 


Schedule  of 
Operations 


The  characteristics  of  all  pesticides  to  be  used  during  the 
rotation  period  are  specified  in  a list.  This  includes 
natural -decay-rate  constant,  solubility,  adsorption  coefficient 
(K  ) , relative  washability  of  the  portion  reaching  plant 
leaves,  and  Arrhenius  constant  if  known  (equation  [117]). 


List  of  Animal  Wastes 

Each  manure  that  is  applied  during  a rotation  cycle  is 
described  in  a parameter  list.  The  list  gives  data  on  percent 
organic  matter,  total  N content,  ammonia  N,  and  phosphorus.  A 
list,  including  nine  common  animal  wastes,  is  incorporated  as 
default  data  in  the  program.  The  user  may  change  any  of  these 
or  add  a tenth,  specifying  the  material's  content  of  organic 
matter,  ammonia,  phosphorus,  and  nitrate.  Inorganic  fertiliz- 
ers are  more  simple;  they  are  specified  by  their  content  of 
nitrate,  ammonia,  and  phosphorus.  The  amounts  of  each  used  are 
given  in  the  operations  schedule,  discussed  below. 


The  lists  described  above  are  for  fixed  characteristics  of  the 
management  operations  that  may  be  repeated  by  several  similar 
operations.  Variable  aspects,  if  any,  of  those  operations  and 
their  dates  are  given  in  the  operations  schedule.  There  are 
separate  schedules  for  tillage  operations,  fertilizations, 
pesticide  applications,  and  irrigation.  A different  schedule 
is  read  for  each  year  of  the  rotation  period,  and  these  opera- 
tions are  repeated  in  the  same  sequence  for  the  number  of 
rotations  required  by  the  simulation  interval. 

The  schedules  needed  include  the  following: 

(a)  Tillages  are  specified  by  number  from  the  tillage  opera- 
tions list  above  and  by  date  of  implementation. 

(b)  Fertilizations  are  specified  by  date  plus  the  rate  of 
application  of  each  of  N,  P,  and  ammonia  (in  lb/acre  or  kg/ha). 
In  addition,  if  animal  wastes  are  the  fertilizer  applied,  the 
number  from  the  above  list  is  given  and  the  type  of  application 
is  specified,  including  depth  of  incorporation  or  injection  if 
used. 

(c)  Pesticides  are  applied  by  specifying  the  pesticide  number 
and  a date,  a code  for  method  of  application,  and  an  applica- 
tion rate  (in  kg/ha  or  lb/acre). 

(d)  Irrigation  schedules  include  specification  of  the  type  and 
amount  of  water  applied  and  the  method  used. 
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Three  types  of  irrigation  may  be  simulated,  as  follows:  (a) 
Sprinkling  on  demand  according  to  soil  water  in  the  root  zone; 
(b)  furrow  irrigation  on  demand  based  on  a soil -water  state; 
and  (c)  ditch  supply  irrigation  for  a specified  period  of  the 
year,  where  water  is  available  on  a rotating  calendrical  basis 
and  with  data  given  on  total  annual  supply.  Furrow  irrigations 
are  routed  along  the  furrow  from  a field  edge  as  described  in 
chapter  5.  For  sprinkler  irrigations,  it  is  assumed  for  sim- 
plicity that  all  applied  water  infiltrates  into  the  soil.  A 
wetted  profile  is  created  and  subsequently  redistributed  just 
as  for  water  entering  the  profile  from  a rainfall. 

Each  fertilization  can  apply  any  or  all  of  N,  P,  and  NH4  in  the 
same  application.  Fertilization  and  cultivation  or  plowing,  or 
any  other  mixture  of  management  operations  can  occur  on  the 
same  date.  However,  each  schedule  must  list  operations  in 
chronological  order. 

Chemigation  is  another  optional  type  of  fertilization  in  which 
fertilizers  or  pesticides  are  applied  with  irrigation  water. 

For  this  operation,  the  date  of  application  of  the  fertilizer 
or  pesticide  must  coincide  with  that  of  the  irrigation,  and  the 
application  must  be  properly  coded  as  a chemigation  method. 
Detailed  instructions  for  specifying  various  management  opera- 
tions is  given  in  the  User's  Manual  (vol . II). 
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9.  TESTING  THE  OPUS  MODEL 


Testing , 
Verification , 
Validation 


Testing  is  commonly  understood  to  be  the  evaluation  of  the 
and  ability  of  a model  to  reproduce  the  measured  characteristics  of 

the  performance  of  an  experimental  prototype.  Given  the  com- 
plexity and  interactivity  of  the  processes  covered  by  the  Opus 
model,  it  is  important  to  understand  not  only  that  no  data  sets 
are  available  to  allow  truly  comprehensive  testing,  but  also 
that  it  is  next  to  impossible  to  know  or  measure  enough  of  what 
is  happening  in  even  a small  agricultural  field,  compared  with 
what  the  model  can  simulate.  Even  with  a simple  model  of  a 
relatively  isolated  process,  the  data  needed  for  rigorous 
testing  should  include  all  the  degrees  of  freedom  of  the  model 
variables  and  input  conditions.  Seldom  are  they  available. 

The  technical  literature  abounds  with  models  of  all  types  that 
have  been  developed  with  appeal  to  generality,  but  are  applied 
to  rather  restricted  or  idealized  conditions. 

The  realistic  appreciation  of  testing  of  a model  such  as  Opus 
requires  a common  understanding  of  the  meaning  of  "testing," 
its  purposes,  and  its  limitations.  First,  testing  of  a model 
should  be  done  only  in  full  understanding  of  the  model's 
intended  use(s).  When  the  Opus  model  is  used  without  calibra- 
tion to  reflect  the  relative  effects  of  alternative  management 
practices,  it  is  improper  to  evaluate  it  on  the  basis  of  the 
accuracy  of  its  uncalibrated  reproduction  of  a given  set  of 
data,  or  to  compare  its  results  with  the  calibrated  results  of 
a different  model. 

There  are  two  general  types  of  testing  in  the  language  of  mod- 
els: verification  and  validation.  Verification  refers  to  the 

process  by  which  it  is  ascertained  that  the  model  is  performing 
as  designed- - that  its  response  is  mathematically  equal  to  the 
equations  on  which  it  has  been  constructed.  This  is  "debug- 
ging." For  a large  complex  model  such  as  Opus,  debugging  is  a 
process  of  diminishing  returns,  and  is  rarely  actually  "com- 
pleted" in  an  economic  time  frame.  Many  years  have  been  spent 
in  the  verification  of  Opus  code. 

Validation  involves  comparison  of  the  performance  of  the  model 
to  that  of  a prototype.  This  aspect  is  the  most  difficult, 
because  it  is  usually  limited  by  the  paucity  (in  extent  and 
intensity)  of  useful  data. 

An  important  point  to  remember  in  testing  any  model  of  a natu- 
ral system  is  the  masking  effect  of  natural  variability  on  the 
accuracy  of  reproduction  of  measured  data.  A perfect  model  of 
some  process  at  a point  on  the  field  may  be  a very  approximate 
model  for  the  combined  action  of  a field  made  up  of  an  assembly 
of  points,  each  having  properties  that  differ  slightly  from 
those  of  its  neighbor.  With  greater  amounts  of  such  spatial 
variability,  the  underlying  point  model  may  be  overwhelmed  by 
the  variation,  if  an  attempt  is  made  to  identify  it  by  the 
measured  performance  of  the  field  as  a whole.  This  does  not, 
of  itself,  negate  the  value  of  the  physically  valid  (point) 
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model,  because  (a)  the  physical  model  is  usually  not  so  masked 
for  every  input  case,  and  (b)  the  problem  indicated  is  not  with 
the  model  but  with  the  complexity  of  the  conditions.  This 
example  represents  many  phenomena  such  as  leaching,  infiltra- 
tion, evapotranspiration , and  surface  erosion. 

The  above  points  are  important  to  remember  when  evaluating  the 
role  of  testing  in  the  determination  of  the  acceptability  of  a 
model.  It  is  just  as  important  to  (a)  understand  the  relative 
merits  of  the  process  components  within  their  disciplines  and 
also  the  testing  that  each  has  already  undergone,  and  (b) 
appreciate  the  structure  of  the  interactions  that  make  up  the 
basis  of  the  model  as  a whole. 


Tests  of  Component  In  looking  at  a model  such  as  Opus,  validation  is  most  appro- 
Processes  priately  applied  on  a component-by-component  basis.  In  this 

regard,  one  should  be  careful  to  fairly  apply  similar  standards 
of  acceptability  to  all  components.  It  is  not  sensible  to 
expect  a high  degree  of  accuracy  in  simulating  the  measured 
distribution  of  sediment-particle  sizes,  for  example,  when  the 
whole  basis  for  the  runoff  phenomena  is  a time  and  spatially 
lumped  method,  which  commonly  exhibits  large  variance  in  pre- 
diction of  daily  amounts.  Interactive  effects  must  also  be 
kept  in  mind  in  validation  efforts.  For  example,  a discrepancy 
that  may  be  observed  between  the  simulated  and  measured  run- 
offs, in  applying  the  model,  may  be  due  to  errors  not  in  the 
traditional  hydrologic  components  (such  as  infiltration  and 
runoff  routing)  but  in  the  simulation  of  the  evapotranspiration 
that  controls  the  initial  conditions. 

The  components  making  up  Opus  each  represent  a considerable 
amount  of  research,  including  validation  of  their  response  to 
the  appropriate  input  effects,  and  it  is  unnecessary  to  ask 
that  such  validation  be  reapplied  for  their  use  in  the  combined 
model.  It  is  proper  to  look  carefully  at  how  the  various 
components  interact,  in  light  of  what  we  think  we  know  about 
how  the  water  cycle  works  on  a small  field.  It  is  also  proper 
to  compare  each  component  with  any  alternate  models  in  each 
respective  discipline  (that  work  with  the  same  level  of  input 
information),  and  to  ask  if  that  component  can  be  improved.  In 
many  cases,  increasing  the  expected  accuracy  involves  an 
increase  in  sophistication,  which  in  turn  requires  more 
detailed  input  information  and  longer  computational  effort. 

This  tradeoff  has  been  considered  in  all  the  components  of  Opus 
and  should  be  kept  in  mind  when  model  testing  is  considered. 


Daily  Weather  Model 

As  indicated  in  chapter  3,  the  WGEN  weather  model  used  in  Opus 
was  extensively  tested  elsewhere  (Richardson  1981) . In  addi- 
tion, several  simulation  runs  of  50  years  with  various  combina- 
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tions  of  options  were  performed,  and  statistics  were  compiled 
to  verify  that  the  model  was  reproducing  the  specified  statis- 
tics . 

Evaporation  from  soil  and  plants  is  calculated  by  modifications 
of  the  method  developed  by  Ritchie  (1972).  The  modifications 
are  minor  and  were  suggested  by  Ritchie  for  EPIC  (Williams  et 
al . 1984).  EPIC  has  been  tested  for  row  crops  (Ritchie  1972, 
Ritchie  et  al . 1976)  and  for  grassland  conditions  (Ritchie  et 
al . 1976,  Pochop  et  al . 1985).  In  addition,  the  Opus  ET  model 
was  applied  to  sites  across  the  United  States  where  pan- 
evaporation data  are  available  for  investigating  corrections 
that  are  applicable  to  account  in  the  average  for  effects  of 
wind  and  relative  humidity  (Opus  User  Manual,  vol . II). 


Rainfall/Runoff  Models 

Extensive  tests  of  the  performance  of  the  daily  hydrology  model 
were  presented  in  the  documentation  to  CREAMS  (Smith  and 
Williams  1980).  As  expected,  because  the  parameter  (the  curve 
number)  that  characterizes  runoff  represents  lumped  results, 
the  performance  of  this  method  is  often  quite  adequate  for  the 
longer  term  although  not  very  good  on  a day-to-day  basis. 

The  erosion  method  accompanying  the  Curve  Number  runoff  method, 
the  MUSLE  model,  was  tested  extensively  elsewhere  (Williams 
1982) . Like  the  Curve  Number  method,  it  should  properly  be 
judged  on  a longer  term  basis;  in  those  applications  it  has 
shown  quite  reasonable  results. 

A simpler  version  of  the  breakpoint  simulation  of  small- 
catchment  hydrology  was  tested  for  use  in  CREAMS  (Smith  and 
Williams  1980) . The  model  in  Opus  makes  several  conceptual 
improvements  over  the  CREAMS  model  without  sacrificing  any  of 
the  basic  methodology  on  which  these  test  results  are  based. 
Further  tests  of  the  analytic/kinematic  infiltration/runoff 
model  were  reported  by  Rovey  et  al . (1977)  and  Smith  (1976, 
1981a) . The  model  for  surface  crusting  effects  used  in  Opus 
was  partly  developed  on  the  results  of  tests  reported  by 
Chevalier  (1984). 


Soil  Water  Flow 

Over  the  last  30  years  or  more,  soil  physicists  such  as  White 
et  al.  (1979)  and  Gilham  et  al . (1979)  have  shown  that 
Richards'  equation  properly  describes  the  one -dimensional  ver- 
tical movement  of  soil  water,  as  long  as  the  general  conditions 
are  appropriate.  More  rigorous  treatment  is  required  for  some 
conditions,  such  as  when  the  escape  of  air  is  severely 
restricted.  Further  testing  of  this  equation  for  use  in  Opus 
should  not  be  necessary.  The  most  common  variation  of  measured 
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data  from  theory  is  the  existence  of  spatial  variation  in  soil 
hydraulic  properties,  rather  than  the  inadequacy  of  the  theory 
or  assumptions  involved  in  Richards'  equation.  Studies  of  the 
effect  of  spatial  variability  (Bresler  and  Dagan  1983)  indicate 
that  model  results  that  ignore  spatial  variability  will  show 
much  less  apparent  dispersion  in  soil-water  content  than  is 
found  in  spatially  variable  soil.  This  bias  should  not  be  a 
serious  detriment,  especially  considering  the  expense  of  data 
collection  and  the  computational  complexity  necessary  to  treat 
spatial  variability  in  the  numerical  simulation.  Results  of 
some  testing  of  the  Separable  Flux  method  of  solution  of 
Richards'  equation  are  presented  in  chapter  4.  Although  the 
increase  in  allowable  time  step  is  a consequence  of  the 
increase  in  layer  size  used,  less  apparent  bias  is  introduced 
than  for  more  conventional  or  linearized  solutions  when  layer 
sizes  are  increased. 


Transport  of  Solutes 

Almost  all  reported  soil -water- transport  data  from  laboratory 
experiments  deal  with  saturated  flow  and  include  diffusion 
processes.  Few  data  are  available  for  unsaturated  flow  from 
intermittent  flows  such  as  those  obtained  in  the  field.  A 
field  data  set  that  allows  limited  testing  comes  from  work  done 
in  Watkinsville , GA  (Smith  et  al . 1978).  This  EPA  report  on 
the  study  describes  the  catchment  and  management  operations  in 
detail.  An  example  of  Opus'  simulation  of  these  data  is  pre- 
sented in  the  followng  pages. 


A relatively  extensive  set  of  test  data  covering  a 3-year  peri- 
od (referred  to  above)  is  provided  by  a study  of  a set  of 
small - field- sized  watersheds  near  Watkinsville,  GA.  This  study 
was  carried  out  by  the  Agricultural  Research  Service  between 
1973  and  1975,  and  forms  the  basis  of  a technical  report  pub- 
lished by  the  Environmental  Protection  Agency  (Smith  et  al . 
1978) . 

The  unique  value  of  these  data  is  the  scope  of  measurements 
made  during  this  period,  including  rainfall;  runoff;  sediment 
production;  soil-water  profiles;  pesticide  and  nutrient  distri- 
butions in  soils;  and  presentation  of  a complete  schedule  of 
tillage  operations,  chemical  applications,  and  resulting  crop 
yields  for  these  fields.  The  experimental  watershed  chosen  for 
this  example  is  called  P-3,  and  a topographic  map  for  this 
field  is  shown  in  figure  42.  The  field  was  planted  with  soy- 
beans each  of  the  3 years  of  the  study,  with  a winter  cover 
crop  of  barley  or  rye.  At  least  one  application  of  two 
pesticides  and  one  application  of  fertilizers  was  used  each 
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year.  A drawback  of  this  data  set  is  the  apparent  lack  of 
consistency  in  the  source  of  raingage  data  for  the  period  of 
the  other  measurements  (Smith  et  al . 1978).  Raingages  were 
replaced  or  data  were  substituted  due  to  mechanical  failure 
during  the  experimental  period.  This  substitution  seems  to  be 
reflected  in  the  simulation  results  shown  below. 


Description  of  Topography 

Figure  43  shows  the  equivalent  furrowed  and  unfurrowed  topogra- 
phy used  in  Opus  to  represent  this  simple  catchment.  The  esti- 
mated mean  slope  length  along  furrows  and  estimated  length  of 
concentrated  flow  path  were  preserved,  along  with  total  field 
area.  The  topographic  map  also  provided  a means  to  estimate 
the  slope  profile  along  the  mean  furrow  flow  path  when  furrows 
control  the  direction  of  runoff  (which  is  most  of  the  time)  as 
well  as  when  furrow  control  is  absent. 

The  measurements  of  outlet  runoff  were  made  with  an  H- flume  in 
the  center  of  an  artificial  dike,  which  forms  a pond.  The 
geometry  of  the  pond  was  taken  from  height  of  the  flume  and 
slope  of  the  land  at  the  outlet.  The  pond- outflow  rating  is 
the  flume  rating,  taken  from  handbooks  (USDA  SEA  1979). 


Description  of  Soil  Horizon 

The  various  soils  mapped  on  this  field  are  similar  in  charac- 
teristics, and  data  to  hydraulically  describe  them  were  found 
from  the  EPA  report  (Smith  et  al . 1978)  and  an  ARS  soil  survey 
(Holtan  et  al.  1968).  When  particular  information  was  unavail- 
able, the  default  options  within  Opus  were  used.  No  extensive 
fitting  or  optimization  was  attempted. 

The  EPA  report  described  the  watershed  soils  as  being  composed 
of  Cecil  sandy  loam  and  sandy  clay  loam,  and  reported 
sand/silt/clay  fractions  for  surface  samples  at  seven  locations 
over  the  catchment  (Smith  et  al . 1978,  table  A3,  p.  173).  From 
this  information,  an  estimated  representative  profile  was  con- 
structed for  use  in  Opus  (shown  in  table  1).  The  upper  and 
lower  horizons  are  sandier  than  the  middle  two,  which  have 
noticeably  higher  clay  contents.  Thus,  a slightly  restrictive 
layer  occurs  at  about  6 inches. 


Management  Actions 

The  schedule  of  tillage,  planting,  harvesting,  pesticide  appli- 
cation, and  fertilizer  application  is  reported  in  appendix  B of 
the  EPA  study  and  is  summarized  here  in  table  2.  The  schedule 
was  translated  directly  into  data  instructions  in  Opus,  and  a 
3-year  simulation  was  performed. 
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Results 


Simulation  of  Runoff 

Table  3 summarizes  the  amounts  of  runoff  and  sediment  produced 
by  the  major  runoff -producing  storms  during  the  3 -year  period 
for  which  data  were  taken.  Like  all  experimental  data,  there 
are  included  some  points  of  quite  doubtful  validity,  such  as 
the  occurrence  of  almost  7 cm  of  rain  without  reported  runoff, 
and  the  reported  runoff  of  almost  6 cm  from  another  7-cm  storm. 
Nevertheless,  comparisons  can  be  made  of  the  comparative  pre- 
diction of  daily  and  breakpoint  models  within  Opus. 

Table  3 contains  only  those  events  with  either  a predicted  or  a 
measured  amount  of  runoff.  Many  storms  are  omitted  in  which 
both  values  are  0.  Thus  the  sample  is  truncated  because  it 
excludes  correct  predictions  in  this  limiting  if  not  trivial 
case.  Often  the  prediction  of  no  runoff  from  storms  is  not 
trivial  but  requires  the  ability  to  use  the  fact  that  the  rain 
intensity  was  very  low.  Indeed,  this  table  shows  many  cases  in 
which  the  Curve  Number  method  predicts  runoff  when  there  was 
none,  and  the  breakpoint  infiltration  model  correctly  predicts 
none  because  of  low- intensity  rain. 

A very  crude  statistical  comparison  of  results  can  be  performed 
using  regression  analysis  on  these  event-runoff  data.  Table  4 
presents  some  of  the  results  obtained  by  this  means.  This 
table  shows  how  differently  each  year's  simulation  performed 
and  also  compares  the  overall  period  results. 

If  one  merely  plots  storm  (or  daily)  runoff  as  a function  of 
rainfall,  it  is  clear  that  a formula  such  as  the  Curve  Number, 
even  considering  initial  soil  conditions,  cannot  explain  the 
variation  in  the  relation.  Indeed,  the  computed  Curve  Numbers 
for  the  storms  on  this  catchment,  based  on  reported  runoff, 
ranged  from  44  to  96.  Another  measure  of  this  is  the  r2  for 
comparison  of  predicted  to  measured  runoff,  using  this  runoff 
model.  Keeping  in  mind  that  the  Opus  daily  model  uses  an 
improved  method  to  relate  changes  in  runoff  to  soil  wetness,  r2 
between  measured  and  CN-predicted  runoff  is  nevertheless  only 
0.49.  Comparative  r2  values  for  all  years  of  breakpoint  simu- 
lation are  0.61.  Forcing  the  intercept  to  0 decreased  this  r2 
to  0.54.  On  a yearly  basis,  the  r2  values  go  up  to  0.77,  0.88, 
and  0.87  for  the  three  successive  years.  Figure  44  shows  the 
breakpoint-model  prediction  for  total  runoff  as  a function  of 
total  rainfall  for  this  data  set,  compared  with  the  reported 
data . 

The  breakpoint  model  simulates  initial  conditions  and  runoff 
hydrographs  for  all  storms  during  a period  when  major  dis- 
turbances affect  the  soil's  capacity  for  intake.  Opus  had 
mixed  success  in  comparison  with  reported  measurements,  and 
some  of  the  simulations  are  shown  in  figures  45  through  53.  As 
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expected,  the  simulation  of  runoff  that  was  a very  small 
percentage  of  total  rainfall  is  the  most  difficult,  and  in 
general,  the  errors  are  largest  in  such  cases.  The  differences 
between  simulated  and  reported  runoffs  are  believed  to  be 
largely  due  to  either  of  two  factors:  quality  of  knowledge 
about  the  true  state  of  the  catchment,  and  quality  of  the 
rainfall - intensity  data. 

Figures  45  through  47  suggest  some  of  the  problems  of  the 
reported  data:  in  these  cases,  apparent  errors  either  in  the 
reported  runoff  or  in  the  reported  rainfall.  It  is  known  from 
comments  in  the  EPA  report  (Smith  et  al . 1978)  that  rainfall 
records  were  not  continuously  good,  and  often  nearby  gage  data 
were  borrowed  to  fill  a record.  Timing  errors  of  a few  minutes 
seem  to  explain  the  cases  in  figures  45  and  46.  Many  hydro- 
graphs were  reported  only  with  a few  data  points,  such  as  in 
figure  45.  Figure  46  suggests  the  loss  of  some  of  the  reces- 
sion data  as  well  as  timing  inaccuracies.  When  the  runoff 
increases  several  minutes  after  the  beginning  of  a rainfall 
pulse,  one  can  often  argue  that  a delay  has  occurred  in  the 
system.  This  should  be  relatively  consistent.  However,  when 
the  runoff  begins  to  decrease  much  earlier  than  the  decrease  in 
the  rainfall  rate,  one  cannot  use  the  same  argument.  This 
seems  to  be  the  case  in  figure  47.  Figure  48  shows  an  excel- 
lent simulation  except  for  a somewhat  clipped  peak  in  the 
reported  hydrograph,  for  which  we  do  not  speculate  a cause. 
Figures  49  and  50  are  common  cases  in  which  runoff  comes  after 
significant  initial  irregular  rainfall  that  prewets  the  soil 
and  is  a severe  test  of  the  somewhat  simplified  soil  physics  in 
Opus.  Nevertheless,  the  results  are  not  severely  in  error 
because  the  runoff  is  only  a small  fraction  of  the  rain  input. 
More  typical  and  successful  simulations  are  shown  in  figures  51 
and  52. 

One  or  two  events  such  as  illustrated  in  figure  53  suggest  that 
the  storage  properties  for  the  catchment  are  larger  than  those 
used  by  the  simulation.  Either  other  processes  are  occurring 
or  the  surface  hydraulic  roughness  is  much  larger.  Overall, 
there  appears  to  be  no  consistent  bias,  and  the  simulation  is 
generally  good.  The  small  but  significant  differences  between 
years  in  statistical  quality  of  the  simulation  and  also  the 
apparent  seasonal  change  in  hydrologic  response  (some  of  which 
is  accounted  for  by  the  crust  model)  suggest  that  the  watershed 
is  very  much  a changing  creature  and  that  the  attempt  to  con- 
tinuously simulate  the  hydrology  remains  a formidable  chal- 
lenge. In  this  respect,  Opus  is  none  too  complex. 


Production  of  Sediment 

With  some  exceptions,  the  results  in  table  3 show  that  most  of 
the  time,  relatively  accurate  simulation  of  runoff  was  accompa- 
nied by  similar  accuracy  in  simulation  of  sediment  production. 
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The  one  major  exception  is  the  storms  on  May  28,  1973.  Given 
the  usual  uncertainty  and  error  in  simulating  a storm's  sedi- 
ment production,  the  relative  success  over  such  a range  of 
storms  is  reassuring  and  adds  to  our  confidence  in  the  model's 
utility  for  transport-simulation  applications. 


Distribution  of  Soil  Water 

The  EPA  report  (Smith  et  al . 1978)  includes  reports  of  soil 
water  averaged  for  depth  increments  and  sampled  at  irregular 
intervals  from  July  1972  to  October  1974.  The  Opus  model 
output  includes  the  option  to  report  soil -water  contents  by 
computational  layer  for  a specified  interval,  and  comparisons 
can  be  made  for  the  dates  of  reported  data. 

In  this  Opus  application,  soil-input  parameters  were  largely 
estimated  by  the  default  option  described  in  chapter  2,  and  the 
only  fitted  parameter  was  the  hydraulic  conductivity  of  the 
surface  horizon.  In  light  of  this,  the  predicted  soil-water 
contents  were  an  excellent  fit  to  the  reported  data,  as  exem- 
plified in  figures  54  to  57.  The  only  significant  bias  appears 
to  be  in  the  porosity  for  the  third  horizon,  where  simulated 
water  contents  are  generally  too  high.  Selecting  a lower  6 
for  that  horizon  could  easily  improve  this  simulation..  Indeed, 
the  use  of  Opus  with  a long  sequence  of  observed  water  contents 
could  be  one  way  to  infer  the  general  shape  of  the  moisture - 
release  curves  of  each  horizon.  One  uncertainty  in  the 
comparison  is  the  actual  time  of  soil  sampling  on  days  when  a 
rainfall  occurred.  This  is  illustrated  in  figure  55,  where 
soil-water  distributions  before  and  after  a storm  are  shown. 

The  relatively  good  fit  at  near- surface  depths  is  also  an 
Indication  of  relatively  good  estimation  of  amounts  and  distri- 
butions of  surface  evaporation  and  root  transpiration  of  water. 


Movement  of  Pesticides 

The  same  EPA  report  (Smith  et  al . 1978)  presents  sample  concen- 
trations, at  various  depths  during  the  study  period,  of  three 
pesticides  applied  at  spring  planting  (see  table  2).  Opus  was 
used  to  study  our  ability  to  simulate  the  measured  concentra- 
tion distributions,  assuming  representative  quality  of  the 
sample  reports.  Both  the  equilibrium  and  kinetic  absorption 
models  were  used,  and  a limited  variation  in  kinetic  rate  and 
interactive  depth  parameters  were  tried. 

The  three  pesticides  were  trifluralin,  paraquat,  and  diphena- 
mid.  Paraquat  is  very  highly  adsorbed  and  moves  only  slightly 
with  water.  The  other  pesticides  are  more  mobile  and  are  a 
better  test  of  the  model ' s ability  to  mimic  actual  soil  trans- 
port. Table  5 lists  the  important  pesticide  parameters  used 
for  the  results  shown  in  the  following  figures. 
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The  EPA  data  (Smith  et  al . 1978)  are  apparently  based  on  one  or 
a few  samples,  although  the  sampling  is  not  discussed.  As 
expected,  this  results  in  considerable  sample  variation,  as  the 
data  show.  Pesticide  concentrations  in  soil  are  reported  for  1 
or  2 years  of  the  3 years  of  hydrologic  data,  so  a limited 
basis  for  comparison  exists.  This  limitation,  coupled  with  the 
year-to-year  variation  in  measured  data  (such  as  apparent  decay 
rates),  diminishes  the  rigor  of  the  tests. 

In  this  test  of  Opus  we  did  not  specify  detailed  vertical  dis- 
tribution of  pesticide  initial  concentrations.  Thus  the 
initial  concentrations  before  first  application  (June  1973)  are 
not  precise  because  the  simulation  starts  in  April  and  some 
mixing  occurs,  but  the  total  amounts  are  small  in  any  case. 

The  situation  after  application  is  much  better,  but  some  small 
limitations  of  the  model  are  apparent.  Opus  (or  any  other 
model)  is  not  able  to  accurately  predict  the  nonuniform  redis- 
tribution of  pesticides  at  a particular  spot  in  the  field  after 
a cultivation  event.  The  mixing  is  not  uniform  either  in  space 
or  with  depth,  as  Opus  assumes.  Neither  is  Opus  able  to 
separately  treat  the  surface  degradation  of  pesticides,  as 
distinguished  from  root-zone  degradation  as  a whole.  Thus 
Opus'  results  are  handicapped  in  predicting  very  precisely  the 
surface  concentrations,  which  are  crucial  to  predicting  runoff 
concentrations.  Moreover,  although  Opus  can  predict  the  con- 
centration of  mobile  material  in  the  surface  (10  mm)  due  to 
soil  evaporation,  it  appears  that  Opus  is  not  able  to  accurate- 
ly predict  the  higher  surface  values  because  of  the  extreme 
gradients  near  the  surface.  Opus  simulates  the  total  amount  in 
the  top  10  mm,  but  probably  much  more  occurs  in  the  top  milli- 
meter than  in  the  tenth  millimeter.  Thus,  Opus  is  expected  to 
underpredict  some  amounts  carried  off  in  sediment. 

Figure  58  compares  simulated  and  reported  concentration  pro- 
files of  diphenamid  at  14,  25,  and  87  days  after  application. 
The  kinetic  adsorption  model  was  used.  The  largest  error  is 
related  to  the  initial  overestimate  of  pesticide  in  the  upper 
100  mm.  The  loss  with  time  over  this  season  is  well  simulated. 
Of  course,  not  all  measured  data  were  simulated  well.  Figure 
59  shows  an  anomalous  period  of  10  days  in  which  the  measured 
amount  of  diphenamid  at  the  surface  fell  far  below  that  simu- 
lated, reversing  a previous  steady  trend.  Other  anomalies  are 
indicated  in  table  6,  which  compares  amounts  of  pesticides  in 
runoff  (both  dissolved  and  attached  to  sediment)  to  reported 
amounts,  for  the  sampled  storms  during  the  period  of  measure- 
ment. An  expected  sampling  variability  is  involved  in  the 
reported  data  here,  because  only  a few  samples  were  taken 
during  the  runoff  and  because  concentrations  will  be  time 
variable.  In  almost  all  cases,  a lower  value  of  runoff 
pesticide  is  associated  with  a lower  value  of  sediment- 
associated  pesticide.  But  diphenamid  on  July  27,  1974,  shows 
much  higher  sediment  pesticide  with  lower  runoff  pesticides. 
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Also,  on  June  11,  1975,  the  reported  soil  concentrations  are 
much  lower  than  those  simulated,  but  the  reported  sediment 
concentrations  of  diphenamid  are  high.  Nevertheless,  the 
results  are  generally  good  in  simulating  the  general  amounts  of 
pesticides  in  the  runoff  and  sediment  several  months  after 
application.  In  most  cases,  simulation  of  the  amount  of 
paraquat  in  sediment  was  as  good  as  the  prediction  of  sediment 
mass . 

In  general,  the  kinetic  model  simulated  the  persistence  of 
pesticides  at  the  surface  better  than  did  the  equilibrium 
adsorption  model.  The  difference  is  typified  by  figure  60. 

The  difference  was  not  large,  which  may  indicate  that  these  two 
mobile  pesticides  are  rapidly  desorbed  and  adsorbed. 

The  equilibrium  and  kinetic  adsorption  models  give  very  similar 
results  for  paraquat  because  this  pesticide  is  so  highly  bound 
to  the  soil.  In  this  example  the  kinetic  model  is  usually  more 
accurate  in  predicting  the  amount  of  pesticide  on  sediment. 
Measured  data  for  dissolved  amounts  of  the  two  more  mobile 
pesticides  give  mixed  results  in  comparisons  of  the  equilibrium 
and  kinetic  models.  Neither  model  is  consistently  superior  in 
simulation. 

Figure  61  shows  samples  of  Opus'  simulation  of  paraquat  distri- 
bution in  the  surface  soil  during  1973  and  1974  in  comparison 
to  reported  measurements.  The  results  from  June  21,  1974, 
illustrate  the  apparent  persistence  of  paraquat  at  the  surface 
compared  to  the  simulated  persistence.  A small  error  in 
simulating  the  actual  mixing  efficiency  for  such  a persistent, 
immobile  pesticide  can  make  all  the  difference  in  the  simulated 
amount  at  the  surface.  The  manipulation  of  mixing  efficiencies 
could  easily  match  the  measured  amount,  but  this  is  not  the 
point  of  our  testing  here. 

A notable  difference  was  observed  between  the  rate  of  decay  of 
diphenamid  during  the  growing  season  and  the  rate  that  occurred 
over  winter.  This  change  could  not  be  simulated  by  Opus,  even 
taking  into  account  the  response  to  water  contents  and  tempera- 
ture as  presented  in  chapter  6.  The  difference  at  the  begin- 
ning of  the  growing  season  in  the  second  year  was  notable  in 
relative  terms  but  negligible  compared  to  the  amount  applied. 
Furthermore,  a striking  difference  was  seen  in  the  rate  of 
decay  for  the  growing  season  of  1973  and  the  rate  for  1974, 
which  remains  outside  the  scope  of  Opus  to  simulate.  This  is 
because  the  difference  cannot  be  explained  in  terms  of  the 
processes  on  which  Opus'  simulation  is  based.  Although  this 
may  result  in  a bias  in  the  simulated  amount  of  pesticide  after 
several  months,  the  error  is  quite  small  compared  to  the  amount 
applied . 
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Most  storms  showed  steady  consistency  for  bias  in  surface 
amounts  and  bias  in  runoff  amounts.  This  consistency  somewhat 
validates  the  lumped  estimation  of  the  interaction  of  surface 
water  and  soil  water  during  storms.  As  the  results  indicate, 
although  sampling  at  larger  depths  may  have  been  limited,  no 
data  indicated  that  the  convective  transport  model  in  Opus  did 
not  reasonably  simulate  the  transport  of  pesticides  with  moving 
water  in  this  soil . 


Summary 

Conclusions  that  are  general  to  the  performance  of  Opus  cannot 
be  drawn  from  one  test  example,  but  some  observations  may  be 
made.  Because  Opus  cannot  simulate  the  spatial  distribution  of 
soil  properties  and  because  the  actual  daily  temperature  or 
radiation  data  for  this  period  were  not  available,  expectations 
about  the  accuracy  of  simulation  of  any  storm  event  should  be 
limited.  The  crop  was  simulated  rather  well  in  terms  of 
reported  yield.  Predicted  runoff  correlated  rather  well  with 
measured  values,  most  hydrographs  were  well  simulated,  and 
predicted  soil  moisture  matched  the  reported  measurements  quite 
well.  The  general  movement  of  pesticides  with  soil  water  and 
with  surface  water  and  sediment  was  reasonably  well  simulated. 
The  accuracy  was  limited  as  much  by  our  ability  to  model  the 
complex  pesticide -degradation  processes  and  by  sampling  uncer- 
tainty, as  by  the  model  of  the  movement  mechanisms. 
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Figure  42. 

Plan  of  watershed  P-3  in  Watkinsville,  GA,  one  of  the  set  of 
watersheds  reported  in  the  EPA  study  of  pesticide  transport 
(Smith  et  al.  1978). 
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Figure  44. 

Comparison  of  Opus  simulation  of  runoff  with  reported  data 
for  watershed  P-3.  Only  days  with  runoff  are  shown,  cover- 
ing the  3-year  period  of  simulation. 


Furrowed  flow  path  geometry 


Unfurrowed  flow  path  geometry 
— 79.9m  1 


Figure  43. 

Geometric  representations  of  watershed  P-3  used  by  Opus  for  test  example. 
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Time  (minutes) 

Figure  45. 

Opus  simulation  of  runoff  event  on  August  29,  1974.  From 
the  rainfall  pattern,  a timing  error  is  indicated,  but  the  few 
recorded  data  points  do  not  adequately  define  the  hydro- 
graph. Square  blocks  in  this  and  subsequent  graphs  are 
recorded  rainfall  pulses. 
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Figure  46. 

Opus  simulation  of  runoff  from  P-3  on  June  27,  1974.  A 
timing  error  is  suggested,  but  the  simulation  is  reasonable 
with  this  caveat.  Graph  of  recorded  data  also  points  to  an 
error  on  the  recession;  the  rapid  drop  at  time  1 140  is 
otherwise  inexplicable. 


Time  (minutes) 


Figure  47. 

Comparison  of  Opus  simulation  with  reported  data  for  the 
event  of  September  9,  1 973.  A problem  in  the  data  is 
indicated  by  the  lack  of  correspondence  between  timing  of 
rainfall  and  runoff  increase  (time  1260-1270  and  1300).  Such 
data  problems  confound  validation  of  any  model. 


Time  (minutes) 

Figure  48. 

Simulation  and  reported  runoff  on  July  27,  1974.  This  case  is 
curious  in  regard  to  reported  data:  the  flattening  of  the 
hydrograph  when  rainfall  is  still  high  is  suspect.  Otherwise 
the  fit  is  excellent. 
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Time  (minutes) 

Figure  49. 

A good  example  of  the  case  where  runoff  is  a tiny  fraction  of 
the  rain,  coming  at  the  end  of  a long  period  of  soil  wetting. 
This  is  the  most  sensitive  test  of  runoff  simulation.  Opus 
underestimates  the  reported  amount,  but  this  kind  of  error  is 
not  surprising. 


Time  (minutes)  *10' 


Figure  50. 

Simulation  of  the  storm  of  March  13,  1975,  exemplifying 
sensitivity  of  prediction  of  soil  redistribution  during  storm 
hiatuses,  here  Opus  methodology  apparently  does  not 
accurately  simulate  sustained  wetness  of  the  surface,  and 
fails  to  match  the  second  (and  following)  small  runoff. 


Simulation  of  runoff  event  of  July  8,  1973. 


Time  (minutes) 


Figure  52. 

Simulation  of  the  storm  runoff  of  July  11,1 975.  Here  again, 
runoff  comes  just  at  the  end  of  the  rainfall  burst,  emphasiz- 
ing any  simulation  error.  Opus'  result  is  very  close  to 
reported  values. 
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Figure  53. 

Simulation  of  the  first  event  of  June  27,  1 974.  This  simula- 
tion indicates  much  more  damping  or  hydraulic  roughness  in 
the  actual  system  than  that  simulated  by  Opus.  This  is  the 
only  such  case  in  this  test. 


Volumetric  Moisture  Content 

Figure  55. 

Measured  and  simulated  water  contents  for  June  28,  1974. 
Day-before  values  indicate  active  surface  drying  during  this 
day.  Some  bias  in  simulated  values  at  intermediate  depths  is 
also  indicated. 


Volumetric  Moisture  Content 

Figure  54. 

Comparison  of  measured  and  simulated  soil-water  distribu- 
tion for  July  13,  1973.  Comparison  of  day-before  and  end-of- 
day  simulations  indicates  little  movement  during  the  day. 
Opus  simulations  extended  to  96-cm  depth,  and  maximum 
root  zone  was  assumed  to  be  about  60  cm. 
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Figure  56. 

Measured  and  simulated  water  contents  for  July  29,  1974. 
This  is  further  into  the  growing  season,  and  water  contents 
are  somewhat  reduced  in  intermediate  depths. 
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Volumetric  Moisture  Content 

Figure  57. 

Measured  and  simulated  water  content  distributions  for 
September  4,  1974.  This  is  near  end  of  growing  season. 

Root  zone  is  deeper  and  water  content  is  somewhat  reduced 
down  to  70  cm,  with  active  drying  during  the  day  at  the 
surface. 
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Figure  59. 

Simulated  and  reported  diphenamid  concentrations  over  an 
1 1-day  interval  in  1974. 
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Figure  58. 

Simulated  and  reported  distributions  of  diphenamid  concen- 
trations at  three  sampling  times  during  the  1973  growing 
season. 
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Comparison  of  simulation  differences  for  kinetic  and  equilib- 
rium adsorption  models  on  a typical  day  in  1973.  Relatively 
rapid  transfer  rates  were  used  in  the  kinematic  model. 
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Figure  61. 

Three  examples  of  simulated  and  reported  distributions  of 
soil  concentration  of  paraquat  during  the  1973  and  1974 
growing  seasons.  Paraquat  is  strongly  adsorbed  and  is  thus 
relatively  immobile. 


Table  1. 

Soil  profile  used  for  watershed  P-3 


Depth 

(nun) 

Sand 

(%) 

Silt 

(%) 

Clay 

(%) 

Saturated 

conductivity 

(nun/min) 

e 

s 

9 

r 

0-152 

64 

19 

17 

0.212 

0.32 

0.08 

152-305 

49 

14 

37 

0.169 

0.35 

0.20 

305-610 

46 

16 

38 

0.132 

0.40 

0.24 

610-1520 

51 

27 

22 

0.306 

0.40 

0.21 
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Table  2. 

Management  operations  for  watershed  Watkinsville  P-3,  1973-75 


Management 

operation 

Date 

Depth 

(mm) 

Efficiency 

(mixing) 

Material/Plant 

Amount 

applied 

(kg/ha) 

Fertilizer 

5/22/73 

Ammonium  nitrate-N 

21 

Superphosphate  - P 

19 

Disk  harrow 

5/5/73 

102 

0.5 

Fertilizer 

6/4/73 

Ammonium  nitrate-N 

25 

Superphosphate  - P 

22 

Rolling 

6/4/73 

51 

0.7 

Cultivator 

Pesticides 

6/15/73 

0 

Trifluralin 

1.12 

0 

Paraquat 

1.53 

0 

Diphenamid 

3.36 

Planting 

6/15/73 

76 

0.7 

Soybeans 

10/5/73 

0 

0.1 

Winter  rye 

Harvest 

11/7/73 

Soybeans 

Planting 

1/14/74 

51 

0.1 

Winter  rye 

Harvest 

5/1/74 

Winter  rye 

Fertilizer 

5/22/74 

Ammonium  nitrate-N 

17 

Superphosphate  - P 

15 

Disk  harrow 

5/22/74 

102 

0.5 

Rolling 

5/28/74 

51 

0.7 

Cultivator 

Pesticides 

5/30/74 

0 

Trifluralin 

1.12 

0 

Paraquat 

1.53 

0 

Diphenamid 

3.36 

Planting 

5/30/74 

76 

0.7 

Soybeans 

Rolling 

6/13/74 

51 

0.7 

Cultivator 

Rolling 

7/5/74 

51 

0.7 

Cultivator 

Harvest 

10/18/74 

Soybeans 

Disk  harrow 

10/22/74 

102 

0.5 

Planting 

10/22/74 

51 

0.3 

Barley 

Harvest 

4/15/75 

Barley 

Fertilizer 

5/8/75 

Superphosphate  - P 

15 

Disk  harrow 

5/13/75 

102 

0.5 

Pesticide 

5/28/75 

0 

Trifluralin 

1.12 

0 

Paraquat 

1.53 

0 

Diphenamid 

3.36 

Planting 

5/28/75 

76 

0.7 

Soybeans 

Sweep 

6/16/75 

51 

0.7 

7/9/75 

51 

0.7 
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Table  3. 
Summary  of 
May  1973  to 

amounts  of 
i September 

runoff  and 
1975 

sediment  for  ' 

W’atkinsville  watershed  P-3, 

Date 

(mm/dd/yy) 

Storm 

rainfall 

(cm) 

Measured 

runoff 

(cm) 

Simulated 

breakpoint 

(cm) 

Runoff 

CN 

(cm) 

Measured 

sediment 

(kg) 

Simulated 

sediment 

(kg) 

5/23/73 

2.21 

0 

0 

0 

10 

0 

5/28/73 

4.83 

1.6 

1.5 

1.77 

1122 

60 

5/28/73 

4.32 

1.9 

1.1 

1.58 

1618 

30 

6/6/73 

3.94 

0.6 

1.4 

0.44 

405 

410 

6/7/73 

2.21 

0 

0.05 

0.09 

0 

20 

6/9/73 

1.19 

0.5 

0 

0 

728 

0 

7/8/73 

6.43 

2.6 

2.4 

1.21 

1184 

560 

7/14/73 

1.9 

0.6 

0.6 

0 

210 

150 

7/17/73 

0.94 

0.2 

0 

0 

66 

0 

9/9/73 

4.45 

1.5 

1.2 

0.23 

58 

180 

9/13/73 

4.12 

0.7 

0.6 

0.18 

25 

120 

12/5/73 

3.86 

0 

0.6 

0.37 

0 

120 

12/20/73 

2.11 

0.3 

0.2 

0.16 

39 

80 

12/26/73 

2.08 

0 

0 

0.05 

0 

0 

12/31/73 

5.33 

0.9 

0.9 

1.78 

78 

250 

1/20/74 

2.34 

0 

0 

0.05 

0 

0 

2/16/74 

1.85 

0 

0 

0.05 

0 

0 

3/21/74 

1.55 

0 

0 

0 

0 

0 

3/29/74 

1.75 

0 

0 

0.01 

0 

0 

4/4/74 

3.55 

0.26 

0.03 

0.36 

*0 

50 

4/13/74 

2.54 

0.2 

0.5 

0.09 

42 

170 

5/23/74 

6.88 

0 

0.75 

1.63 

0 

30 

6/27/74 

5.33 

0.8 

0.9 

1.59 

196 

240 

6/27/74 

3.3 

1.1 

1.2 

0.99 

557 

450 

7/24/74 

1.17 

0 

0 

0 

0 

0 

7/27/74 

7.7 

1.7 

2.5 

1.53 

537 

490 

8/16/74 

4.45 

0.8 

1 

0.32 

227 

250 

8/29/74 

2.03 

0.2 

0.5 

0 

42 

150 

9/1/74 

1.27 

0.05 

0.04 

0 

16 

20 

11/30/74 

3.56 

0 

0 

0.24 

0 

0 

12/15/74 

3.12 

0 

0 

0.20 

0 

0 

12/19/74 

2.16 

0 

0 

0.05 

0 

0 

1/10/75 

2.59 

0.05 

0.01 

0.08 

22 

10 

1/12/75 

3.12 

0.7 

0 

0.44 

75 

0 

2/16/75 

2.62 

0.1 

0 

0.41 

3.5 

0 

2/18/75 

4.54 

0.5 

1 

1.03 

47 

20 

2/24/75 

2.41 

0.5 

0.001 

0.08 

97 

10 

3/13/75 

10.01 

3.6 

1 

3.23 

258 

250 

3/14/75 

2.54 

0 

0 

0.28 

0 

0 

3/15/75 

3.73 

0 

0.45 

0.82 

0 

1450 

3/24/75 

2.64 

0.3 

0.01 

0.12 

24 

10 

3/30/75 

1.52 

0 

0 

0 

0 

6 

4/2/75 

6.99 

5.8 

2.5 

2.32 

305 

650 

4/14/75 

2.29 

0 

0 

0.02 

0 

0 

5/3/75 

3.94 

0 

0.1 

0.42 

0 

70 

5/7/75 

2.79 

0 

0.03 

0.1 

0 

20 

5/14/75 

1.9 

0 

0 

0 

0 

0 

5/31/75 

3.56 

0 

0.3 

0.06 

0 

200 

6/11/75 

2.67 

0.3 

0.7 

0.63 

73 

220 

6/11/75 

2.54 

0.6 

0.6 

0.59 

168 

160 

9/6/75 

1.79 

0 

0 

0 

0 

0 

9/12/75 

1.85 

0 

0.13 

0 

0 

30 

9/17/75 

3.4 

0 

0.5 

0.12 

0 

60 

9/22/75 

1.55 

0 

0.04 

0 

0 

20 

9/23/75 

4.88 

0.1 

0 

1.25 

3.5 

0 
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Table  4. 

Regression  analysis  of  simulated  runoff 


Sim.  = a+b*Meas . 

Sim.  - 

b*Meas 

a b 

b 

Model 

Sample 

(Intercept)  (Slope)  R2 

Slope 

R2 

Curve 

1973 

0.064 

0.605 

0.49 

0.646 

0.49 

number 

1974 

0.198 

0.799 

0.43 

0.964 

0.37 

1975 

0.256 

0.488 

0.69 

0.554 

0.6 

All 

0.215 

0.521 

0.57 

0.609 

0.49 

Break- 

1973 

0.106 

0.785 

0.77 

0.854 

0.76 

point 

1974 

0.055 

1.269 

0.88 

1.32 

0.88 

1975 

0.12 

0.367 

0.87 

0.398 

0.72 

All 

0.195 

0.503 

0.61 

0.582 

0.54 

Note:  Sim.  = simulated  value, 

Meas . = measured  value . 


Table  5. 


Parameters 

for  pesticides  used  at  watershed  P- 

3,  Watkinsville , GA 

Pesticide 

Basic 

half-life 

(days) 

K 

oc 

(liter/kg) 

Solubility 

(mg/liter) 

t (2> 

z 

Kinetic 
rate  factor, 
(liter/min) 

Diphenamid 

6.9 

400 

260 

0.03 

0.009 

Paraquat 

462.0 

1.6E7 

5.0E5 

0.03 

0.009 

Trifluralin  17.3 

600 

1 

0.03 

0.009 

(1)  Carbon- dependent 

adsorption 

ratio.  See  ch. 

. 4. 

(2)  Surface-water  extraction  ratio.  See  ch.  6. 

(3)  u , defined  for  equation  [36]. 
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Table  6 . 

Summary  of  pesticides  in  runoff  and  sediment  for  major  storms 


In  runoff  (g/ha):  On  sediment  (g/ha) : 

Date  Pesticide  Reported  Kinetic  Equilib.  Reported  Kinetic  Equilib. 


July 

Diphenamid 

16 

3.4 

5.1 

0.75 

0.93 

0.15 

8/73 

Paraquat 

0 

0.025 

0.025 

67 

28 

28 

Trif luralin 

2.2 

1.5 

2.3 

0.11 

0.43 

0.098 

July 

Diphenamid 

2.5 

1.5 

2.7 

0.094 

0.067 

0.1 

14/73 

Paraquat 

0 

0.0043 

0.0043 

11 

6.5 

6.5 

Trif luralin 

0.28 

0.76 

1.4 

0.015 

0.049 

0.079 

June 

Diphenamid 

3.4 

1.5 

2.8 

0.23 

0.17 

0.11 

27/74 

Paraquat 

0 

0.009 

0.009 

31 

13.2 

13.2 

Trifluralin 

- 

0.66 

1.52 

- 

0.09 

0.085 

July 

Diphenamid 

0.0017 

0.051 

0.11 

0.088 

0.0027 

0.00045 

27/74 

Paraquat 

0 

0.009 

0.0091 

18 

1.4 

1.4 

Trifluralin 

- 

0.08 

0.22 

- 

0.0044 

0.0013 

May 

Diphenamid 

- 

4.0 

9.9 

- 

2.0 

1.7 

31/75 

Paraquat 

- 

0.0016 

0.0016 

- 

11 

11 

Trifluralin 

- 

0.9 

2.2 

- 

0.53 

0.55 

June 

Diphenamid 

1.1 

3.2 

7.2 

0.062 

0.5 

0.46 

11/75 

Paraquat 

0 

0.0061 

0.0061 

4.5 

15.2 

15.2 

Trifluralin 

- 

1.2 

2.8 

- 

0.22 

0.26 

(-)  indicates  data  not  measured  or  reported. 
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